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Abstract 


This final report presents new algorithms for target detection, feature extraction, 
and image formation with the synthetic aperture radar (SAR) technology. For target 
detection, we consider target detection with SAR and coherent subtraction. We also 
study how the image false alarm rates are related to the target template false alarm 
rates when target templates are used for target detection. For feature extraction from 
SAR images, we present a computationally efficient eigenstructure-based 2D-M0DE 
algorithm for two-dimensional frequency estimation. For SAR image formation, we 
present a robust 
of radar targets 


parametric data model for estimating high resolution range signatures 
and for forming high resolution SAR images. 
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1. Introduction 


This final report presents new algorithms for target detection, feature extraction, 
and image formation with the synthetic aperture radar (SAR) technology. 

In Chapter 2, we consider target detection with SAR and coherent subtraction. 
We shall show with some limited experimental data that the coherent subtraction 
technique may be used to suppress outliers and obtain approximate Gaussian distri- 
butions for clutter and noise. We shall also derive generalized likelihood ratio (GLR) 
detection algorithms that may be used with SAR images that have Gaussian distribu- 
tions. We shall analytically compare the performance of a) a single pixel detector, b) 
a detector using complete knowledge of the target signature information and known 
orientation information, c) a detector using incomplete knowledge of the target sig- 
nature information and known orientation information, d) a detector using unknown 
target signature information and known orientation information, and e) a detector 
using unknown target signature information and unknown orientation information. 

In Chapter 3, we study how the image false alarm rates are related to the target 
template false alarm rates when target templates are used for target detection. In 
particular, we shall show a simple way of determining the probability of false alarm 
of a target template when a low constant false alarm rate is desired for an image and 
when the image size is much larger than the size of the target template. 

In Chapter 4, we present a computationally efficient eigenstructure-based 20- 
MODE algorithm for two-dimensional frequency estimation or feature extraction from 
SAR images. We derive the theoretical performance of the 2D-M0DE estimator 
and show that it is asymptotically statistically efficient under either the assumption 
that the number of temporal snapshots is large or the signal-to- noise ratio is high. 
Numerical examples showing the performance of this algorithm and comparing it with 
the computationally efficient subspace rotation algorithms are also given. We show 
that the statistical performance of the 2D-M0DE algorithm is better than that of 
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the subspace rotation methods, whereas the amount of computations required by the 
former is usually no more than a few times of that needed by the latter. 

In Chapter 5, we present a robust parametric data model for estimating high res- 
olution range signatures of radar targets and for forming high resolution SAR images. 
This paper also presents an estimation algorithm for the data model. The algorithm is 
referred to as the APES (Amplitude and Phase Estimation of a Sinusoid in unknown 
colored noise) algorithm. We shall describe how the APES algorithm can be used to 
estimate range signatures and to form SAR images. We shall show, with both numer- 
ical and experimental examples, that our modeling and estimation approach yields 
better resolution and lower sidelobes than the conventional nonparametric FFT (fast 
Fourier transform) method. We shall also show that our approach is more robust 
than modeling the radar data as a certain number of complex sinusoids in noise and 
estimating the frequencies, amplitudes, and phases of the sinusoids with one of the 
best sinusoidal parameter estimation methods. 

Each of the afore-mentioned chapters is self-contained with its own introductions, 
formulations of the problems of interests, approaches, conclusions, and references. 

The results we present in this report are obtained with the partial support from 
WL/AARA, Wright Laboratory, Wright Patterson Air Force Base as a subcontract 
from the Ohio Aerospace Institute. Our work is also supported in part by the 1993 
and 1994 AFOSR Summer Faculty Research Programs, by the National Science Foun- 
dation, and by the Goran Gustafsson Foundation. 

Those who contributed to this report include Dr. Jian Li, Mr. Syed M. Rahman, 
Dr. Petre Stoica, Mr. Edmund G. Zelnio, and Mr. Dunmin Zheng. 
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2. Target Detection with Synthetic Aperture Radar and Coherent 

Subtraction 


2.1 Introduction 

Synthetic aperture radar (SAR) technology may be used to detect radar targets 
of interest. High resolution SAR technology is especially useful for detecting small 
radar targets embedded in strong ground clutter such as in foliage. In this paper, 
we shall consider target detection algorithms that may be used with high resolution 

SAR. 

Target detection from SAR or optical images has been considered by many au- 
thors. For example, Novak, Burl, and Irving [1] considered target detection with a 
polarimetric SAR. The three output images of the polarimetric SAR are first processed 
by a polarimetric whitening filter, which is derived by assuming a K-distribution for 
clutter and noise. The output image of the filter is next used with a two-parameter 
detector for target detection. The target detection in [1] is performed one pixel at a 
time even though the target may occupy more than one pixel; i.e., even though the 
target size may be larger than the resolution of the SAR image. 

Reed and Yu [2] considered generalized likelihood ratio target detection from a 
sequence of optical images, which are first preprocessed by removing local means so 
that the clutter and noise will approximately have the Gaussian distribution. In [2], 
each target in an image is described by a completely known template or signature 
with an unknown gain, which is a scalar. The algorithm, however, may not work well 
with SAR images. For example, for a target in foliage, the SAR target signature may 
change due to the interaction between target and surrounding clutter. 

Stotts [3] considered detecting several dim targets in an image simultaneously. 
The image is also first preprocessed by removing local means so that the clutter and 
noise will approximately have the Gaussian distribution. Each dim target is described 
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by a known template or signature with an unknown gain, which is a scalar. Stotts 
has shown that simultaneous detection of multiple targets may perform better than 
separate detection of individual target. In this paper, we shall extend this idea of the 
simultaneous detection of multiple targets to the detection of a target with multiple 
pixels in a SAR image. 

We shall derive generalized likelihood ratio (GLR) detection algorithms that may 
be used with multiple SAR images that are obtained with coherent subtraction or 
have Gaussian distributions. We shall show with some limited and well calibrated 
experimental data that we may eliminate Gaussian outliers of the clutter and noise 
through coherent subtraction between two complex SAR images of the same area of 
interest. One of the two SAR images is assumed to be target free and the other is to 
be sought for the presence of target. In each image obtained with coherent subtrac- 
tion, the target of interest is modeled with a target template, which is large enough 
to cover the entire target. The size of the target template and the number of complex 
unknowns in it are determined by the knowledge of the target orientation informa- 
tion and the amount of target signature information known to the detector. Using 
this unifying framework, we shall analytically compare the performance of a) a single 
pixel detector, b) a detector using complete knowledge of the target signature infor- 
mation and known orientation information, c) a detector using incomplete knowledge 
of the target signature information and known orientation information, d) a detector 
using unknown target signature information and known orientation information, and 
e) a detector using unknown target signature information and unknown orientation 
information. We shall derive the probability of detection and the probability of false 
alarm of each detector. To achieve a constant false alarm rate (CFAR), each detector 
threshold is simply a function of the dimensional parameters of the detection problem 

and the desired probability of false alarm. 

In Section II. we discuss the effects of Gaussian assumption on target detection. In 



Section III, we formulate the target detection problem. In Section IV and Appendices 
A, B, and C, we assume that the statistics of the clutter and noise are known and 
derive a GLR detector under the assumption and discuss its performance. In Section 
V and Appendices D and E, we assume that the statistics of the clutter and noise are 
unknown and present a more practical GLR detector. In Section VI, we apply the 
more practical detector to the experimental data. Finally, Section VII contains our 
conclusions. 

2.2 Gaussian Assumption 

A key problem in radar target detection is the description of the statistical prop- 
erties of radar clutter and noise. In general, radar clutter and noise do not satisfy 
the conditions of a Gaussian distribution. Many statistical models, such as the well- 
known log-normal, Weibull, and K-distributions, have been proposed to describe the 
clutter and noise statistics. Although these distributions may provide better statis- 
tical models for the clutter and noise than the Gaussian distribution, the detectors 
that are derived based on these models may be very complicated and may involve 
an expensive multidimensional search over the parameter space. It is also difficult 
to analyze the performance of these detectors. It is difficult, if not impossible, to 
derive the probability of detection and probability of false alarm for such a detector. 
It is thus difficult to analytically determine the correct threshold for a prescribed 
probability of false alarm. 

Alternatively, the SAR images may be preprocessed so that the Gaussian assump- 
tion for clutter and noise is approximately valid. For example, Hunt and Cannon 
[4] and Reed and Yu [2] considered removing local means as such a preprocessing 
technique. The detectors obtained from Gaussian assumptions may avoid the multi- 
dimensional search over the parameter space. It is also quite tractable to analyze the 
performance of such detectors and many Gaussian assumption based results exist in 
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the literature. Yet if the distribution of the preprocessed images is still non-Gaussian, 
then the detectors obtained under the Gaussian assumption will suffer performance 
degradation. The selection of a detector, therefore, must balance these tradeoffs. 

In this paper, we consider target detection by assuming the radar clutter and noise 
in a S AR image have a Gaussian distribution. Such a SAR image may be obtained by 
taking the difference between two SAR images of the same area, with one assumed to 
be target free and the other is to be sought for the presence of target. We show below 
that the coherent subtraction method can suppress Gaussian outliers when the SAR 
images are obtained by moving the radar antenna along a fixed rail. These rail SAR 
images can be used in many applications including monitoring the environmental 
changes of a cetain area of interest. For an airborne SAR, however, we do not know 
how the coherent subtraction method may perform. The presence of errors such as 
the flight path errors and the existence of atmospheric turbulence can degrade the 
performance of the coherent subtraction. 

The limited and well calibrated experimental data we have are obtained by ERIM 
(Environmental Research Institute of Michigan) with a portable rail SAR that has 
foliage penetration capabilities [5]. More specifically, a horizontal 36 feet long alu- 
minum truss was used to support a rail and an antenna carriage. The data was 
collected by moving an antenna along a fixed rail to obtain the synthetic aperture. 
The rail SAR is an FM-CW radar system based on an HP-8501 network analyzer 
[o]. Our results are, therefore, tentative and the effectiveness of coherent subtraction 
needs to be further studied with more extensive experimental data obtained under 
more operational conditions. 

The data we shall use were obtained when both the transmitter and receiver 
of the SAR are horizontally linearly polarized. Figures 2.1(a) and (b) show the 3- 
dimensional (3-D) plots of the magnitudes of two complex SAR images obtained with 
two identical synthetic apertures. The frequency band used for the image is between 
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400 and 1300 MHz and the depression angle is 30°. We have downsampled the original 
images presented in [5] by a factor of two in range and by a factor of six in cross range 
since the original images are oversampled. The range and cross range resolutions in 
the SAR images are 0.34 and 3 meters, respectively. Figure 2.1(a) shows the 3-D plot 
with foliage only. The peaks in the figure correspond to the radar returns from tree 
trunks. Figure 2.1(b) shows the 3-D plot of a target in foliage. The target is a pickup 
truck rotated 24° counterclockwise from end-on. Figure 2.1(c) shows the 3-D plot 
of the magnitude of the coherent subtraction between the two complex SAR images. 
Figure 2.2 is similar to Figure 2.1 except that the truck is broadside. Thus the target 
return in Figure 2.2 is much stronger than in Figure 2.1. 

We note that coherent subtraction can effectively suppress the large clutter returns 
due to tree trunks, which result in false alarms and causes CFAR detectors to fail. 
The large returns left in Figures 2.1(c) and 2.2(c) are due to the target and its 
surroundings. For example, the darkened peak to the left of the target in Figure 
2.2(c) occurs in Figure 2.2(a) but not in Figure 2.2(b). This result occurs because the 
main response of a tree comes when the radar energy bounces off the tree onto the 
ground and returns to the radar or vice versa. The presence of the target interrupts 
this path and results in no tree return in the image in Figure 2.2(b). Subtracting 
the image in Figure 2.2(a), which contains the tree return, from the image in Figure 
2.2(b), which does not contain the tree return, results in a tree return not canceled 
by the subtraction process thus yielding the darkened peak. The presence of the 
darkened peak in Figure 2.2(c) is additional information showing the presence of a 
target because its presence is due to the interaction between the target and clutter. 
This information may be especially useful when the target return is weak. 

Figure 2.3 shows the magnitude of the correlation coefficient of the clutter and 
noise (i.e., the target-free) pixels in Figure 2.1(c) as a function of the spatial distance 
(as measured by the number of pixels) between two pixels. An unbiased autocorrela- 
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tion estimator [6] was used to estimate the correlation coefficient. We note that the 
clutter and noise pixels are approximately independent of each other. 

Consider next the test of normality of the clutter and noise pixels before and after 
coherent subtraction. Let x„ i = 1, 2, • • • , Q, denote a 2 x 1 vector containing the 
real and imaginary parts of the ith pixel used for the test. Let 


Q 


where (-) r denotes the transpose, 


and 


- (x, - x) T S 1 (xi-x) , 

(2.1) 

1 Q 

(2.2) 


(2.3) 


For true Gaussian random vectors, the variance of # is approximately equal to 4 / (Q - 
l) 2 . Thus the test of normality of x, may be performed by comparing the sample 
variance of # with 4 f(Q - l) 2 [7]; i.e ;j by comparing how close c is to 1, where e 
denotes the ratio between the sample variance of /?, and 4 /(Q - l) 2 - For the clutter 
and noise pixels in Figure 2.1, the e is equal to 7.97 before coherent subtraction. This 
large e is caused by the large tree trunk returns. After coherent subtraction, however, 
the e is equal 1.10 and is much closer to 1. The histogram of the clutter and noise 
pixels after coherent subtraction can also be shown to match a zero-mean complex 

Gaussian probability density function well. 

Finally, since broadcast stations often occupy low frequency bands, SAR images 
of the same area may be formed with two or more separate frequency bands instead 
of one wide frequency band. The separate bands may be chosen to avoid the jamming 
from the radio stations. The clutter and noise pixels in different images obtained with 
different frequency bands can also be show n to be approximately independent of each 

other. 
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The problem formulation below will take these results into account. 

2.3 Formulation of the Target Detection Problem 

Consider J high resolution SAR images obtained with coherent subtraction in 
which a target may be present. The J images may be obtained from polarimetric 
SAR, different frequency bands, and/or different aspect angles. For each image, the 
target may be modeled with a template consisting of Nj pixels, j = 1,2, • ■ • , J. The 
shape of the templates may be arbitrary and the templates may consist of areas that 
are not connected. Among the Nj pixels, Kj ( I<j < Nj ) pixels are assumed to be 
deterministic and arbitrary unknown complex scalars that correspond to the bright 
returns of target scatterers. The remaining Nj - Kj pixels are assumed to contain 
clutter and noise only; i.e., they correspond to the areas of the target that does not 
generate radar returns. 

The locations of the I<j pixels may be assumed known, not completely known, 
or unknown. If the locations of the unknown scalars are not known exactly, the 
dimension Kj (and Nj correspondingly since Nj > Kj) may be increased to include 
all possible locations of the scalars. If the locations are unknown and we only know 
the approximate size of a target, we may choose a template large enough to cover the 
target and assume that all pixels in the template are unknown, i.e, Kj = Nj. Yet 
increasing Kj decreases the detection performance due to the increased number of 
unknowns in the target template. This result will be shown in the following sections. 

Let z denote an N x 1 vector consisting of the pixels of such templates in the 
presence of clutter and noise, where N = N\ + N2 + ■ • ■ + Nj. Under hypothesis H\, 
the target presence hypothesis, the z may be written 

z = Sb + n. (2.4) 

The b is the I\ x 1 vector consisting of the I\ deterministic unknown complex scalars 
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in the target templates, where K = K\ + Ki + ■ • • + Kj- The b may be written 

b= [bf b T 2 bj] , (2-5) 

where bj, j = 1,2 , • • • , J, is the I<j x 1 vector consisting of the Kj deterministic 
unknown complex scalars in the target template of the j th image, or the j th template. 

The S is a full-rank NxK matrix describing the locations of the unknown scalars. 
Only one element in each row and each column of S is one and the remaining of the 
elements are zero. Thus we have 

S h S= 1 k , (2.6) 


where denotes the complex conjugate transpose and I/c denotes the identity 

matrix of dimension K. The S may be written 



where Sj, j = 1, 2, • • • , J, is a full-rank Nj x I<j matrix describing the locations of 
the unknown scalars of the j^th target template and 



( 2 . 8 ) 


The n denotes the N x 1 clutter and noise random vector and may be written 



(2.9) 


where nj are Nj x 1 clutter and noise vectors of the j th image and are assumed 
zero-mean complex Gaussian with covariance matrix ctj and are statistically in- 
dependent of each other. Under hypothesis Hq , the target absence hypothesis, the z 
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may be written 


z = n. (2.10) 

The problem of interest is to develop a CFAR detector for the data model (2.4). We 
shall consider the generalized likelihood ratio (GLR) target detection algorithms for 
the purpose. We shall consider both a detector, where the clutter and noise variance 
is assumed known, and a more practical detector, where the clutter and noise variance 
is unknown. The effects of the dimensional parameters such as N, A', and J on the 
performance of the detectors will also be considered through performance analysis of 
the detectors. 

We remark that the above detection problem is a generalized version of the ap- 
proach of considering one pixel at a time and the approach of using the complete 
knowledge of the target signature. When we consider detection by using one pixel 
at a time, we have Nj = Kj = 1, j = 1,2, ■■■,«/. When the complete knowledge 
of the target signature is known except for a complex gain, the S j becomes the 
signature vector s j whose elements are arbitrary except that its Euclidean norm is 
constrained to be 1. Also, the b j becomes the unknown complex scalar gain bj and 
Kj = 1, j = 1, 2, * - • , J. The elements of s j describe the locations and relative return 
strengths and phases of the scatterers of a target. 

We next remark on the knowledge of S j and s j, which may be obtained through 
experiments or simulations with software packages similar to XPATCH [8]. The target 
signature vector s j may be easily altered by the surrounding environment of a target, 
but the location matrix Sj is more robust to the target surroundings. For a given 
target, a set of signature vectors or location matrices may be obtained for different 
orientation angles of the target relative to the radar. These vectors or matrices may 
be used with GLR detectors to form a filter bank. 

For example, consider detecting the target shown in Figures 2.1 and 2.2. Since 
J = 1, for example, we drop the subscript j for simplicity. For this example, we may 
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use four different S’s to form a filter bank: (i) when the target is oriented broadside 
(or +180° due to symmetry), (ii) when the target is oriented end-on, (iii) when the 
target is oriented head-on, and (iv) for the rest of the target orientations. We note 
from Figure 2.1 that for Case (iv), the target return consists of only one bright pixel 
and hence we should let S = 1, which results in a single pixel detector. 

When a filter bank is used for target detection, the probability of false alarm 
may be increased as compared with using a single filter; i.e., a single Sj or S ; , for 
target detection. When a set of S^’s is used, for example, the increase will depend 
on how similar the Sj’s are. The more non-overlapping elements in S j s, the larger 
the increase of the probability of false alarm. \et the increase is no more than the 
number of filters times the probability of false alarm of using a single filter. It appears 
that the exact probability of false alarm of using a filter bank for target detection 
cannot be determined analytically in general and must be obtained with Monte-Carlo 
simulations. In other words, although a constant false alarm rate can be achieved by 
using a filter bank, we may not know exactly what the achieved false alarm rate is. 

The design of a filter bank is more flexible using S j than using s j. We may increase 
the template size Nj and the number of unknown template pixels Kj to reduce the 
number of filters in the filter bank. Increasing the number unknown template pixels 
Kj has the effect of decreasing the sensitivity of the target detector to target signature 
information. Increasing the template size Nj has the effect of decreasing the sensitivity 
of the target detector to target aspect. By selecting appropriate Kj and Nj for each 
image based on the amount of the a priori knowledge about the target, we may use 
a single filter instead of a bank of filters for target detection. The next paragraph 
illustrates these issues with several examples using different combinations of Kj and 
Aj- 

Assume that we have J = 1 image and we drop the subscript j for simplicity. 
Consider, for example, the four scenarios shown in Figure 2.4. In Figure 2.4(a), 
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we assume we have the complete knowledge of the target signature except for the 
unknown gain. Also, the target orientation is assumed known. For this case, we have 
J<[ = ]_ unknown in the target template and we should use §6 instead of Sb to model 
the target return. In Figure 2.4(b), we assume that we have the incomplete knowledge 
of the target signature; i.e., we assume that we know the locations of the unknown 
scalars that represent the bright returns from the target scatterers. The number of 
the unknown scalars in the target template is assumed to be A =10 for this case. 
We also assume that the target orientation is known and the number of pixels in the 
target template is, say, N = 100. For this case, then, S is a 100 x 10 matrix and b 
is a 10 x 1 vector. In Figure 2.4(c), we assume that the target orientation and the 
target approximate size are known. For this case, all pixels in the target template are 
assumed unknown and we assume that there are K = 100 unknowns in the target 
template. For this case, then, S is a 100 x 100 identity matrix and b is a 100 x 1 
vector. In Figure 2.4(d), we assume that only the approximate size of the target is 
known. Since the target orientation is unknown, we choose a target template that is 
large enough to cover all possible target orientations and all pixels in the template 
are assumed unknown. We assume that there are K = 1000 unknowns in the target 
template. For this case, then, S is a 1000 x 1000 identity matrix and b is a 1000 x 1 
vector. In the following sections, we shall show how the target detection performance 
is affected by the use of these target templates. 

2.4 A CFAR Detector Based on Known Clutter and Noise Variance and 
Its Performance 

We shall present below a generalized likelihood ratio (GLR) detector under the 
assumption that the clutter and noise variance in each image; i.e., a), is known, and 
also present its performance. This detector is referred to as Detector A below and its 
performance is easy to analyze. 
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( 2 . 11 ) 


It is shown in Appendix A that Detector A has the form 

J v Hc.cH„. Hi 

y z i j >* 

U ff i 

The threshold parameter 7 is determined according to a given probability of false 
alarm and thus the detector is a CFAR detector. 

It is shown in Appendix B that the probability of false alarm of Detector A is 


K — 1 K—l—k 

Pf=E TT? i — exp(— 7 ). 




( 2 . 12 ) 


We note that the Pf depends only on /v, the sum of the numbers of unknown param 
eters in the target templates, and the threshold parameter 7 . The Pf is independent 
of Nj, the sizes of the templates. For a given probability of false alarm, the 7 in 
the detector (2.11) is obtained with (2.12). It is also shown in Appendix B that the 
probability of detection of Detector A is 


00 k 6 ' t+ ^A 1 

= f- 7 )ga g (i + ff-i-*)! 


i+K—l—k 


(2.13) 


where 

3 h^b 

« = E7r- < 2 - 14 ) 

i=i 

The 5 is the sum of the signal-to-clutter-and-noise ratios (SCNRs) of the templates. 
Note that 8 is the sum of the signal-to-clutter-and-noise ratios of the non-zero pixels 
in all templates. We also note that the Pq is also independent of the template sizes 
Nj. The P D depends on K, SCNR 8, and the threshold 7 . 

In most of the following examples, we consider Pf = 10 -10 for the target templates. 
This is because we intend to achieve a false alarm rate of 10 4 per km 2 , where km 
denotes the kilometer. (ARPA’s goal is to achieve a false alarm rate of 10 ~ 3 per 
km 2 .) For a SAR image with range and cross range resolutions of 0.34 and 3 meters, 
respectively, the probability of false alarm of a single pixel detector is approximately 


15 



lO -10 . When the target template size is much smaller than the image covering the 
one km 2 area, the template false alarm rate should also be approximately 10 10 [9]. 

To compare the performance of using different target templates with the Detector 
A, let us consider the effect of I< on the performance of Detector A. Figure 2.5 shows 
the probability of detection as a function of SCNR for different A when P F - 10 10 . 
When we have one image; he., J = 1, for example, the four performance curves shown 
in Figure 2.5 could correspond to the four scenarios shown in Figure 2.4. It is shown 
that for the given P F and a fixed SCNR <5, the P D of the detector in (2.11) decreases 
as IC increases. We note from Figure 2.5 that the best performance occurs when we 
use the complete knowledge of the target signature except for the unknown gain and 
when the target orientation is known. The worst performance occurs for the case of 
unknown target signature and unknown target orientation. This result occurs due 
to the non-coherent integration loss [10, p.71] since as may be seen from (2.43) in 
Appendix B, the energies of the template pixels are summed up non- coherently. 

Note also that for SCNR = 15 dB in Figure 2.5, the average SCNR per pixel is 
15 dB for I< = 1 and is 5 dB for I< = 10. Thus Figure 2.5 also shows the effect of 
smearing a given amount target energy among more pixels and hence lowering the 

SCNR per pixel. 

Figure 2.6 shows the extra SCNR needed to achieve P D = 0.5 for different proba- 
bilities of false alarm P F . We note that we have similar curves for different probabil- 
ities of false alarm P F . The extra SCNR needed decreases slightly as P F decreases. 
We also note that the extra SCNR needed increases slowly as A increases. 

Although it is the best to use the complete knowledge of the target signature for 
target detection, the target signature may not be completely known for SAR images 
and even when known, the signature may change due to the interaction between the 
target and clutter and other factors. The change of signature may result in severe 
detector performance degradation. Let Sj and s_, be the assumed and true target 
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signature vectors, respectively, where both s j and Sj have Euclidean norm 1. It is 
shown in Appendix C that the probability of false alarm for this case is the same as 
(2.12) with I< = J. The probability of detection for this case has the form of (2.13) 


with K — J and 




(2-15) 


1 

J 


where Pj (0 < p 3 < 1) are the SCNR loss factors as a result of the signature mismatch; 



Figure 2.7 shows the SCNR loss as a function of pi when J = 1. 

We now make the comparison between using the incomplete knowledge of the 
target signature information and the approach of using the complete knowledge of the 
target signature information except for the unknown gain. The comparison may be 
made most easily with J = 1. For J = 1, we drop the subscript of p and I\. Consider 
the example where J = 1 and we use an incomplete target signature described by I\ = 
20 unknown parameters. The extra SCNR needed to achieve the same probability of 
detection as when using the complete target signature information ( I< = 1) is about 3 
dB, as shown in Figure 2.6. If we know that the mismatch between the assumed and 
true target signatures is small (possibly through experiments), or more specifically, 
if the mismatch will not result in p < 0.5, then it is better to use the complete 
knowledge of the target signature information in the detector. Otherwise, it is better 
to use the incomplete signature information and assume K = 20 unknowns in the 
target template for this example. 

The comparison between considering one target template at a time and one pixel 
at a time is also clear. As shown in Appendix C, for the same target SCNR, the 
performance of assuming the complete knowledge of the target signature is the same 
as ft. — \ j = 1.2. , J, in our problem formulation. For J = 1 and a template 

with I\ - 20 unknown parameters, for example, the extra SCNR needed to achieve 
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the same probability of detection as for K = 1 is about 3 dB, as shown in Figure 
2.6. Then for the case where the target return is a single bright pixel, the net loss of 
using the target template is about 3 dB, as compared with using one pixel at a time. 
If the template SCNR 8 is at least 3 dB larger than the SCNR of the highest pixel 
in the template, then, using the target template model is better than the single pixel 
detector. For the best possible case where all K = 20 pixels in the target template 
have equal magnitude, the template SCNR is about 13 dB more than the individual 
pixel SCNR. Thus the net gain of using the target template for target detection is 
about 10 dB because of the 3 dB loss due to the increased number of unknowns. 

We now compare the performance of using a bank of target templates at a time 
and using a single pixel at a time for target detection. Note that the maximum 
probability of false alarm of using a filter bank for target detection is the number of 
target templates or filters times the probability of false alarm of using a single target 
template or filter. Thus we may set the probability of false alarm of each target 
template to be the desired probability false alarm divided by the number of filters. 
The probability of false alarm of the single pixel detector, however, is the same as the 
the desired probability false alarm. Thus we encounter an additional SCNR loss when 
using a bank of target templates for target detection due to the decreased probability 
of false alarm for each target template. Yet when the number of filters in a filter bank 
is small, this SCNR loss is negligible. When the number of filters is 4, for example, 
this SCNR loss is much less than 1 dB. (See Figure 2.10 for example.) 

Finally, let us consider the advantages of using multiple SAR images for target 
detection. Consider the example where J = 2 and AT = AT = 20. For the best case 
where the target template SCNRs for both images are assumed the same, the net 
gain of using both images for target detection is about 2 dB because of the 1 dB loss 
due to doubling the number of unknowns, as shown in Figure 2.6. 
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2.5 A CFAR Detector Based on Unknown Clutter and Noise Variance 
and Its Performance 


In the previous section, we have studied the performance of a CFAR detector that 
assumes that the clutter and noise variance in each image is known. In practice, 
however, the clutter and noise variances are unknown. We present below a detector 
for this practical situation. This CFAR detector is referred to as Detector B. Our 
approach is similar to the one developed by Kelly [11] for target detection with a 

phased array airborne surveillance radar. 

The Detector B we shall present utilizes both primary and secondary data of a 
SAR image for target detection. The data vector z,-, from which the target presence 
is sought, is referred to as the primary data. For the jth SAR image, the secondary 
data vectors are denoted Zj (l), ■ ■ -, z ,-(£;)• These vectors are assumed to be 

target free; i.e., they represent the target free background of the jth SAR image. 
They are assumed to have the same statistics as Zj, the subvector of the primary 
data vector z, under hypothesis H 0 and are statistically independent of each other 
and z. The secondary data are useful for estimating the clutter and noise variance in 


Detector B. 

It is shown in Appendix D that Detector B has the form 

,Hr. , V-A rrH(" 1 (Lj+l)7Vj 


J 

n 

J=1 


zfej + Egi gfXOMO 


Hi 

H 0 


(2.17) 


[zf (Tv, - S]sfJ 7j + e5i (oi(0 J 

The threshold parameter £ is determined according to a given probability of false 

alarm and thus Detector B is a CFAR detector. 

In general, we simply set the number of the secondary data vectors to be equal to 

each other; i.e., we let 

L l =L, = ~- = L J ±L (2.18) 
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Taking the natural logarithm of the ( Lq + l)st root of (2.17) yields 


r] = ^ln 

3 - 1 


zjz, + Tti £ (0^(0 


_zf (l Wj - S,S f) Zj + Zt, zfVMO. 


Hi 

<7, 

H 0 


(2.19) 


where 7 = (ln£)/[(Z 0 + l)Hj]. 

It is shown in Appendix E that except for <7 = 1, there are no closed form ex- 
pressions for the probability of detection Pd and probability of false alarm P F ■ It is 
shown that under hypothesis Hq, 77 has the pdf 


fMH 0 ) = f m {v\Ho) * fmin |Ho) * • ' • * fnjiv |H 0 ), 

where * denotes the convolution and for j = 1, 2, • • • , J, 

. . n > [ eX V(v) Z 1 l Aj ~ 1 iMl + K J Z *) ! 

Jvpvi oj - ^ _ 1)! exp[(Mj + Kj - 1)77] ’ 


( 2 . 20 ) 


with 


( 2 . 21 ) 


( 2 . 22 ) 


Mj = Nj{L 0 + l)-Kj. 

Under hypothesis iZi, 77 has the pdf 

/„(77|tfi) = f m (v |Hi) * f n {nWi) *■■■* fiAv IHi), (2.23) 

where for j = 1, 2, • • • , J, 

, , ^ SfeMv) - l] A Z+»- 1 (M j + Kj + i - 1)! . 

fM H ') - ex P(-^) L i!(M . _ 1)!(A ^. + j _ 1 )! exp[(M, + Kj + z - 1)77]’ 1 ' J 


where 6, is the SCNR of the jth template. 


* 

Cj — 2 


(2.25) 


The P F and P D may be obtained numerically from f v (v\Ho) and fn (v\ R i)i respec- 
tively; i.e., 

Pf = J fr,{rj\H 0 )dT}, (2.26) 
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and 


(2.27) 


r co 

Pd= Mv\Hi)dr). 

The fvinlHo) given in (2.20) and f n (rj \H X ) given in (2.23) may be obtained more 
efficiently by first calculating the characteristic functions of fr,j(r)\H 0 ) and f V} {T)\H x ) 
with FFT (Fast Fourier Transform). For a given probability of false alarm; i.e., to 
achieve CFAR, the 7 in the detector (2.19) is obtained with (2.26). 

We remark that the above analysis also holds when L 0 is not an integer but a 
rational number such that NjL 0 , j = 1,2, ■■■,./, are integers since the number of 
target free pixels in each image does not have to be the multiples of the size of the 
corresponding target template. The detector (2.17) or (2.19) may be changed slightly 
to accommodate the fact that Lo is not an integer. 

Consider next J = 1, where we have the closed form expressions for P D and P F . 
For this case, we drop the subscripts of z, z (/), 8, M, K, and N. Taking the ( L + l)st 
root of (2.17) yields the detector for J = 1: 

+ £L z g (/)z(Q _ >V (2.28) 

Z" (Iv - SS* ) Z + Ef=i "0 

It is shown in Appendix E that the probability of false alarm of the above detector is 


K - 1 


Pf=Z 

k = 0 


(C - 1)*-*- 1 

£ M+K-k-l 


M + K-k-2 
I< - k - 1 


(2.29) 


We note that for J = 1, the P F of Detector B depends on the dimensional parameters 
M and K and the threshold parameter (. For a given probability of false alarm, the 

£ in the detector (2.17) is obtained with (2.29). 

It is also shown in Appendix E that the probability of detection of the above 

detector when J = 1 is 


A/-1 


P D = 1 - exp (-6/C) Yj 

k—0 


M-k - 1 


E 


' M + I< - k - 2 
^ K + i- 1 


(g/CV(C-i) A ' +i 

t ! SM+K-k - 1 ' 


(2.30) 
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We note that for J = 1, Pd depends on the dimensional parameters M and A, the 
SCNR <5, and the threshold (• 

We now examine how the different dimensional parameters affect the performance 
of the detector in (2.17). These effects of the parameters may be most easily explained 
with the case of J = 1. The effects are then generalized to the case of J > 1. Let 
j — 1 and consider M = NL + N — A', the total number of target free pixels in 
the primary and secondary data vectors. The larger the number of target free pixels 
M, the better the estimate of the clutter and noise variance, and hence the closer 
the performance of Detector B to that of Detector A. This result may be observed 
from Figure 2.8, which shows the probability of detection as a function of SCNR for 

different M when J = 1, K = 2 and Pf = 10 

Consider next the effect of A', the number of unknown parameters in the target 
template, on the performance of Detector B. We first explain that the larger the K, 
the more number of target free pixels M is needed by Detector B to achieve similar 
performance as Detector A. We shall consider J = 1 and drop the subscript j for 
convenience. As shown in Appendix B, we may rewrite Detector A in (2.11) as 

= > 7 '. (2-31) 

A H 0 

Under hypothesis H 0 , z A has the complex Gaussian distribution with zero-mean and 
covariance matrix I*. Under hypothesis H u z A has the complex Gaussian distribution 
with mean b/cr and covariance matrix I*. Thus under hypothesis H 0 , the mean and 
variance of r?' are 1 and 1/A', respectively. Under hypothesis H u the mean and 
variance of rj' are 1 + 8/K and 1/A' + 8/K, respectively, where 8 is the SCNR given 
in (2.14). As shown in Appendix E, we may rewrite Detector B in (2.17) as 

„ gjWH >7" (2.32) 

771 " [zgzs + £f=i z"(/)z(/)] /M & 

The statistical properties of the numerator of r/" is the same as the properties of rj’ in 
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(2.31). Since under both hypotheses, the z (/) has the complex Gaussian distribution 
with zero-mean and covariance matrix I N , the mean and variance of the denominator 
of 7 )[ are 1 and 1/M, respectively. We note that for large K and small M, the 
performance of (2.32) will be affected by the variance 1/M of the denominator of 77" 
Thus for large K, M must also be large in order for Detector B in (2.32) to achieve 
similar performance as Detector A in (2.31) that assumes that the clutter and noise 
variance is known. We found through numerical examples that for \ <K< 1000, 
the performance differences between Detectors A and B are similar for different A 
when M is proportional to K 2/ 3 . Figure 2.9 shows the probability of detection as a 
function of SCNR for different K when J = 1, M = [4SA 2/3 J and P F = lO” 10 , where 
(X] denotes the integer part of x. The figure also shows the performance of Detector 
A for comparison. Note that Detectors A and B have similar performances. 

Figure 2.10 shows the probability of detection as a function of probability of 
false alarm; i.e., the receiver operating characteristic of Detector B, for different I< 
and SCNR when J = 1, M = US, and P F = IO -10 . We note that the 5 dB 
change in SCNR or the change in the number of unknowns I< from 1 to 10 has a 
significant effect on the P F for a given P D . This significant change is due to the 
exponentially decreasing tail associated with the assumed Gaussian model for the 
clutter distribution. 

The discussions for the case of a single image ( J = 1) may be easily generalized to 
the case of multiple images ( J > 1)- Consider J = 2, for example. To achieve similar 
performance as Detector A, the number of target free pixels in each image must be 
large for Detector B. Figure 2.11 shows the probability of detection as a function of 
SCNR for different Mi = M 2 when J = 2, I\\ = I< 2 = 2 and P F = 10 10 . We note 
that the larger the Mi and M 2 , the closer the performance of Detector B to that of 

Detector A, which is expected. 
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2.6 Target Detection with Experimental Data 

We consider first the performance of Detector B when used with the experimental 
data shown in Figure 2.1. Figure 2.12 shows the detection results before and after 
coherent subtraction when M = 48, N = I< = 1, and P F = lO' 10 . Figure 2.12 shows 
the generalized likelihood ratio (GLR) obtained with the left side of (2.28) when it 
is above the detection threshold and zero when it is below the threshold. We note 
that before coherent subtraction, although the presence of the target is detected, the 
large tree trunk returns also result in a false alarm. Thus in the presence of large tree 
trunk returns, Detector B is no longer a CFAR detector and its probability of false 
alarm is also determined by the number of large tree trunk returns. With coherent 
subtraction, however, the false alarm due to the large tree trunk return is eliminated 
and Detector B is truly CFAR. 

Consider next the experimental data shown in Figure 2.2. We note that the target 
occupies more than one pixel. Let us assume that the target orientation is known 
and use an incomplete target signature described by I< = 20 unknowns. (Note that 
the incomplete target signature may be determined by simulation softwares such as 
XPATCH [8] or by experiments in a laboratory. They should not be determined from 
images from which the target presence is sought. When the target SCNR is small, 
determining target templates from such images is impossible.) Then the template 
SCNR for the target is approximately 35 dB for the data shown in Figure 2.2(c). The 
largest pixel SCNR for the target is about 25 dB. Thus compared with using each 
single pixel for target detection, using the template with K = 20 unknowns results 
in a net gain of about 7 dB because of the 3 dB loss due to the increased number of 

unknowns. 

Finally, we could add noise to the data shown in Figure 2.1 to simulate a weak 
target in clutter and noise. Assume that with the added noise, the template SCNR 
for the target is approximately 20 dB and the largest pixel SCNR for the target is 
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about 10 dB. Then from Figure 2.9, we note that the probability of detection of using 
the target template is 1 while the probability of detection of using the single pixel 
detector is approximately 0. 

2.7 Conclusions 

We have considered target detection with synthetic aperture radar. We have de- 
rived generalized likelihood ratio (GLR) detection algorithms that may be used with 
SAR images that are obtained with coherent subtraction or have Gaussian distri- 
butions. Through performance analysis, we have analytically compared the perfor- 
mance of a) a single pixel detector, b) a detector using a complete knowledge of the 
target signature information and known orientation information, c) a detector using 
an incomplete knowledge of the target signature information and known orientation 
information, d) a detector using unknown target signature information and known 
orientation information, and e) a detector using unknown target signature information 
and unknown orientation information. 
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Appendix A: Derivation of Detector A 


Under hypothesis , the probability density function (pdf) of the complex Gaus- 
sian random vector z may be written 


J 1 

/i( 2 ) = n N . 2N- ex p 


j = 1 


(Zj - S j 6 i ) g (z J - Sj-fy) 


(2.33) 


Under hypothesis H 0 , the pdf of the complex Gaussian random vector z may be 
written 

j i r '7^'z-i 

(2.34) 


J 1 

/o(z) = n -v ex p 


i=i 


TT^O* 2 ^ 


zf Zj 


j j 


max b /i 


The generalized likelihood ratio has the form - . Maximizing /i with respect 

Jo 

to bj yields 

b, = Sf Zi, (2.35) 


and 


J 1 

max /i(z) = n^exp 




7r N J<j 2N J 


zf (Iat,. ~ SjSf ) Zj 


(2.36) 


where we have used (2.8). Then the generalized likelihood ratio test becomes 


n«p 


’•zfSjSfzj 


i=i 




(2.37) 


Taking the natural logarithm of both sides of the above test yields the optimal detector 

J -HccH-. h x 

E J r - " £ 5 7, (2-38) 

U °i 

where 7 = ln£. 


Appendix B: Performance of Detector A 

Consider the GLR test in (2.11). Let Zj = z j/aj. Then under hypothesis H 0 , 
Zj has the complex Gaussian distribution with zero-mean and covariance matrix Iv,. 
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Under hypothesis Hi, Zj has the complex Gaussian distribution with mean Sjbj/crj 
and covariance matrix Tv, . The GLR test in (2.11) may now be written 


E z f SjSf z, > 7 . 

j= 1 


w 0 


(2.39) 


Since SfSj = I,v j5 there exists a unitary matrix Uy such that 


VjSj = 


0 


(2.40) 


Let z j = UjZj. Then under hypothesis H 0 , z, j still has the complex Gaussian dis- 
tribution with zero-mean and covariance matrix I,y, • Under hypothesis H\, has 
the complex Gaussian distribution with mean \J jSjhj /aj and covariance matrix Tv,- 
The GLR test in (2.39) may now be written 

X^SfUjSjSfUfS,- > 7 . (2.41) 

4 — 1 


Let 



(2.42) 


where I j A and z jB are Kj x 1 and (Nj — Kj) x 1 vectors, respectively. By using (2.40) 
and (2.42), (2.41) may be written 


V = 



(2.43) 


Under hypothesis H 0 , ZjA has the complex Gaussian distribution with zero-mean 
and covariance matrix Tv,- Under hypothesis Hi, Zja has the complex Gaussian 
distribution with mean bj/<7j and covariance matrix Ia' 7 - Thus under hypothesis Ho, 
2 r) has the central x 2 distribution with 2 K degrees of freedom since K = AT + AT + 


h AT; i.e., 


fMH 0 ) 


77 A ' l exp(-ii) 

(K - 1)! • 


(2.44) 
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Then the probability of false alarm is 


/*CO 

Pf= /*(»? \H 0 )dr) 

J ' y 


K - 1 A'-l-fc 


(2.45) 


Under hypothesis Hi , 2 tj has the noncentral x 2 distribution with 2 A degrees of free- 
dom and noncentrality parameter 6 = J2j = i b^bj/a 2 ; i.e., 


») = ex P(- { “ ,) g7!(; + /f--Tjr 


(2.46) 


Then the probability of detection is 

TOO 

Pd= fv(v\ H i) dr l 

oo t+A'— 1 ^i+A'— 1 — k 

= E (j + 

Appendix C: Effect of Target Signature Mismatch 


(2.47) 


When the target signatures Sj are used in the GLR test in (2.11), (2.11) becomes 


J '7 ^ O 

U °j 


(2.48) 


where s^ij = 1. Under hypothesis H 0 , z j has the complex Gaussian distribution with 
zero-mean and covariance matrix <7 2 Tv, • Under hypothesis H \ , Zj has the complex 
Gaussian distribution with mean s jbj and covariance matrix 

Through similar transformations as in Appendix A, the test (2.48) may be written 


v = £M 2 jV 


(2.49) 


Under hypothesis H 0 , z jA has the complex Gaussian distribution with zero-mean and 
variance 1. Under hypothesis H u z jA has the complex Gaussian distribution with 
mean sfsjbj/<Tj and variance 1. Thus under hypothesis Hq, 2 t] has the central x 
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distribution with 2 degrees of freedom. Under hypothesis H\, 2r] has the noncentral 
X 2 distribution with 2 degrees of freedom and noncentrality parameter 


We note that in the absence 
in (2.50) becomes 


* = E hf- 


j=i a i 


of target signature mismatch; i.e., when 


« = E 


M 


9 ? 

i=i <U 


(2.50) 
s j, the 8 

(2.51) 


which is the SCNR for the case where the target signature is known completely. 


Appendix D: Derivation of Detector B 


The derivations below are both a specialization (due to the uncorrelated clutter 
and noise pixels) and a generalization (due to multiple images) of [11]. Under hy- 
pothesis H u the probability density function (pdf) of the complex Gaussian random 
vectors z j, z,(l), • ■ •, z ;(!;), j = 1, 2 , • ■ ■ , J, may be written 


/i ( z j) z j(i)» ” ■ > 


i,2,-w)=n 

i=x 


N 

7 


: exp 



(2.52) 


where 


Tii* 


Xj + 1 


G 


( Zj - S,6 J ) H (z i - S fr) + £ z f (0 z i(0 


;=i 


(2.53) 


Under hypothesis if 0 , the pdf of the complex Gaussian random vectors Zj, Zj(l), • • 
z j = 1, 2, • • • , J, may be written 


fo ( z j> z j(l), ■ • • ,Zj(Lj),j = 1,2, • ■ ■ , J) — I] 

;=i 


" 7 r^ 7 exp 

71 ^ CT; J 



Lj +1 


(2.54) 


where 


T 0] = 


1 

Lj + 1 



(2.55) 
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max a 2 ... b /i . . . 

The generalized likelihood ratio has the form *’ ’ J ’ • Maximizing f 0 with 

max„!... 5 2 /o 


respect to a] yields 


Thus 


«T? = —T 0j . 


3 Ni 
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(L } +l)Nj 
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Maximizing f\ with respect to a 2 yields 
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= —T,y 


Thus 
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Minimizing Tu with respect to bj gives 

b i = Sfzy 


and 


minTiy = 


Lj + 1 


zf (i*. - S,Sf) zj + £zf (0^(0 

1=1 


Then the generalized likelihood ratio test becomes 

zfzj + Efei z f(0 z i(0 


J 

n 

j=i 
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Appendix E: Performance of Detector B 


Hi 

It- 

HO 


(2.56) 


(2.57) 


(2.58) 


(2.59) 


(2.60) 


(2.61) 


(2.62) 


The derivations below are again both a specialization (due to the uncorrelated 
clutter and noise pixels) and a generalization (due to multiple images) of [11]. 

J=L L 

Consider first J = 1 and the GLR test in (2.28). Let 




(2.63) 


30 



Note that t}[ is a monotonically increasing function of r)'. Then (2.17) may be rewritten 

z H SS H z >£_i. (2.64) 


Vi = 


z H (Ijv - SS* ) z + £f=i "o 


Let 


and 


Uz 


z — 


2(0 = 


Uz(0 


l = 1,2, • ■ ■ , L, 


(2.65) 


( 2 . 66 ) 


where U is a unitary matrix that satisfies (2.40). Then under hypothesis Ho, z 
has the complex Gaussian distribution with zero-mean and covariance matrix Ijv. 
Under hypothesis H u z has the complex Gaussian distribution with mean USb/cr 
and covariance matrix I N . Under both hypotheses, the z (/) has the complex Gaussian 
distribution with zero-mean and covariance matrix Ijv- By using (2.65) and (2.66), 
(2.64) becomes 

,;= g „ g ffl- 4 -fc-1, (2-67) 

where 


+ ZL 2"(0^(0 4. 


z = 


zb 


( 2 . 68 ) 


with z A and z B denoting K x 1 and (N-I<) x 1 vectors, respectively. Under hypothesis 
Ho, z a h as the complex Gaussian distribution with zero-mean and covariance matrix 
I A -. Under hypothesis H u z A has the complex Gaussian distribution with mean b/a 
and covariance matrix I K -. Under both hypotheses, z B has the complex Gaussian 
distribution with zero-mean and covariance matrix I n-k- Let 


t = z%z A , 


(2.69) 


and 


T = zfz B +^Z H (/)z(/). 


(2.70) 


;=i 
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Then (2.67) becomes 


(2.71) 


t <{ C - l) T - 
^0 

Under hypothesis H 0 , 2 1 has the central x 2 distribution with 2 1< degrees of freedom. 
Under hypothesis Hi, 2 1 has the noncentral x 2 distribution with 2 A degrees of free- 
dom and noncentrality parameter 8 = b H b jo 2 . Under both hypotheses, 2 r has the 
central x 2 distribution with 2(iVL + N - A') = 2M degrees of freedom. The < and r 
are independent of each other. Thus under hypothesis H 0 , = Mtf (I< t) has 

the central F-distribution with 2 K and 2 M degrees of freedom. Under hypothesis 
Hi, rj[ M j K — Mtf(Kr ) has the noncentral F-distribution with 2A and 2 M degrees 

of freedom and noncentrality parameter 8. 

The probability of false alarm may now be calculated as 


P F = P{t > « - l)r\H 0 ) 

r oo r roc 

= \ f,(t\Ho)Jt 
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fr{r)di 
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oo I /■ oo 

«-l)7 

(M + K - fc - 2)! 


- X ^w+A'-fc-i (/<• _ k - 1 )!(M - 1)! ’ 


(2.72) 


where / t (*|tf,) and f T \ Hi (r) denote the conditional probability density functions of t 
and t, respectively, under hypothesis Hi, i = 1,2. 

The probability of detection may be calculated as 

P D = 1- [°°\[ fr{r)dT f t {t\Hi)dt 

Jo |A/(c-i) J 

= i-«pM)L (M-i-ijr 
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X i! 
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Using the finite sum expression derived by Kelly in [12], we have 


a/-i M—k—i ( m + K - k - 2 
Pd = 1 - exp (-8/0 Y X! 
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J>h 

Consider below J > 1 and the GLR test in (2.19). Let (2.19) be rewritten as 


* = I>(gi + 1) < ln ^’ 

i= i K ° 


where 


zf s j s fzj 


(2.75) 


(2.76) 


qi zf (In, - S,Sf ) z, + zf m(0 ' 

Note that the qj is similar to t]\ in (2.64). Thus under hypothesis Ho, qjMj/hj has 
the central F-distribution with 2 1\ 3 and 2 Mj degrees of freedom; i.e., 

f ( | in gf-^ + ^-l)! 7) 

J V mHo) - _ !)!(! + qj )M ]+ K i ■ 

Under hypothesis H \ , qjAIj/Kj has the noncentral F-distribution with 2 Kj and 2 Afj 
degrees of freedom and noncentrality parameter 5j\ i.e., 

fMth) = expM,) g i!(M /l ? l) ! ( ir /^^| , (i + +?j )Mi^+." (2 ' 78) 

Note also that <71, 92? ■ • <?j are statistically independent of each other. Let 

(2.79) 




i n (?j + 1 )> i — 1? 2, • ■ • , J- 


It is easy to show that 


(( 1 rr v [exp (jh) - l] Aj 1 (Mj + Kj - 1)! 

fniWWo) - ^ K . _ _ i)l exp[(Mj + Kj - 1 )%]’ 


(2.80) 


and 


5i[expfe) - + *)+■- 1)! n 

•Mwlffi) - ex P( s i) g ,i( iW , _ 1)!(/C, + i - l)!exp[(A/, + Kj + i - 1)%] ' 

The 771, 772, • • •, rjj are also statistically independent of each other. Thus 

/„0jI« 0 = /„(-?!%) * /*(>?l«>) * • ■ • * Mim, (2-82) 
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and 


fMHi) = fm (v\Hi) * fmtn \ H 0 * • • • * fnj(v\Hi)- 


(2.83) 
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Figure 2.2: 3-dimensional plots of the magnitudes of complex SAR images when the 
target is oriented broadside and frequency band between 400 and 1300 MHz. (a) 
Foliage only, (b) Target in foliage, (c) After coherent subtraction. 
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Figure 2.3: Magnitude of clutter and noise correlation coefficient as a function of 
spatial delays measured by the number of pixels. 


k=i k = io 



Figure 2.4: Target template examples for (a) complete knowledge of target signature 
information except for the unknown complex gain and known orientation informa- 
tion, (b) incomplete knowledge of target signature information and known orienta- 
tion information, (c) unknown target signature information and known orientation 
information, and (d) unknown target signature information and unknown orientation 
information. 
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Figure 2.5: Probability of detection of Detector A vs. signal-to-clutter-and-noise ratio 
for different number of unknown parameters K with P F = 10 -1 °. 



Figure 2.6: As compared with K = 1, extra SCNR needed for Detector A to achieve 
P D — 0.5 vs. K for different probabilities of false alarm P F . 
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Figure 2.7: Signal-to-clutter-and-noise ratio (SCNR) loss for Detector A as a result 
of target signature mismatch vs. SCNR loss factor when J = 1. 



Figure 2.8: Probability of detection of Detector B vs. signal-to-clutter-and noise ratio 
for different M when J = 1, I\ = 2, and Pf = 10 10 . 
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Figure 9 9- Probability of detection vs. signal-to-clutter-and-noise ratio for different 
K when J = 1 , M = L48/i 2 / 3 J, and P F = lO" 10 . The dashed and solid lines are for 

Detectors A and B, respectively. 



Figure 2.10: Receiver Operating Characteristic of Detector B for different h^a.nd 
different signal-to-clutter-and-noise ratio when J = 1, M = 118, and P F = 10 
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Figure 2.11: Probability of detection of Detector B vs. signal to-clutter-and-noise 
ratio for different M\ = M 2 when J = 2, L\ = Z -2 = -Mi/2 = M 2 /2, Ai = A 2 - 2 
and = 10“ 10 - 


Figure 2.12: Target detection (end-on) with J — A — A — 1, M 
P F = 10“ 10 . (a) Before coherent subtraction, (b) After coherent subtraction. 



3. On Image and Template False Alarm Rates When Using Target 

Templates for Target Detection 


3.1 Introduction 

Due to the availability of high resolution synthetic aperture radar (SAR) and 
optical images, ground targets such as trucks and tanks are often described by target 
templates with each template consisting of one or more pixels in an image. As a 
result, target templates are used for target detection [1, 2]. It has been shown in [1] 
that when a target is described by more than one bright pixel, using an appropriate 
target template for target detection may give better performance than using one pixel 
at a time. 

For a given probability of false alarm of a target template, generalized likelihood 
ratio detectors may be derived to achieve the constant false alarm rate (CFAR) [1, 2]. 
Yet achieving CFAR for each target template may not be what a radar designer is 
asked to achieve. For example, a radar designer may be asked to achieve a constant 
false alarm rate per kilometer squared. For this case, the size or dimension of an image 
that covers an area of one kilometer squared is determined by the image resolution. 
Thus a CFAR detector should be derived to achieve CFAR for the image rather than 
for the target template. 

However, there does not appear to exist an analytical expression that shows how 
the false alarm rate of a target template is related to the false alarm rate of an image 
except for a couple of special cases. One of the special cases is when the size of the 
image is equal to the size of the target template. For this case, the two false alarm 
rates are the same. The other special case is when the target template is a single 
pixel and the pixel false alarm rate is small. For this case, the image false alarm rate 
is approximately equal to the template or pixel false alarm rate multiplied by the 
number of pixels in the image. 
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In this letter, we shall show with computer simulations how the false alarm rate 
of a target template is related to the false alarm rate of an image when we wish to 
achieve an image CFAR, which is much less than 1. We shall show that for the case 
where the radar clutter and noise pixels are independently and identically distributed 
and the image size is much larger than the target template size, the image false alarm 
rate is also approximately equal to the template false alarm rate multiplied by the 
number of pixels in the image. This case is often encountered in practice. For the 
example where each image covers an area of one kilometer squared, the image size 
is easily much larger than the target template size for most ground targets such as 
trucks and tanks. 

3.2 Simulations Results 

The simulation results below are obtained for the following target detection prob- 
lem [1]. Let z denote aniVxl vector consisting of the pixels of a target template in the 
presence of clutter and noise. Under hypothesis H u the target presence hypothesis, 
the z may be written 

z = Sb + n. 

The b is the I< x 1 vector consisting of K deterministic unknown complex scalars 
in the target templates. The S is a full-rank N x I< matrix describing the locations 
of the unknown scalars. Only one element in each row and each column of S is one 
and the remaining elements are zero. The n denotes the N x 1 clutter and noise 
random vector and is assumed zero-mean complex Gaussian with covariance matrix 
a 2 I, where I denotes the identity matrix. Under hypothesis H 0 , the target absence 
hypothesis, the z may be written 

z = n. 

It has been shown in [1] that the generalized likelihood ratio detector for the above 
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target detection problem has the form 

z H SS H z l 1 

2 <?• 

a 1 Ho 

The threshold parameter 7 is determined according to a given probabdity of false 
alarm of the target template Pp(Target Template), where 

K - 1 K—\—k 

P F ( Target Template) = Y (ft _ 1 _ ky. exp ( _7 ' - 

To detect the presence of a target from a SAR or optical image, the target tem- 
plate, to which 3.2 is applied, is shifted both horizontally and vertically and one pixel 
at a time across the entire image. Note that as the target template moves across the 
entire image, the areas covered by the template may overlap. When more than one 
detections occur over a group of overlapping areas of the image, all detections are 
combined together as one detection since targets do not overlap. 

In the simulation examples below, the images from which the presence of a target is 
sought are assumed to be L x L square images. The target templates are also assumed 
to be square templates with N = K. Both the image and the target template false 
alarm rates are assumed to be much less than 1. Also, the template false alarm rate 
is assumed to be so small that the probability of two or more false alarms occuring 
over non-overlapping areas of an image is negligible as compared to the probability 
of one false alarm per image. Thus in the simulation examples below, the image false 
alarm rates are computed as the number of images that contain false alarms divided 
by the total number of independent images used in the Monte-Carlo simulations. 

Figure 3.1 shows the image false alarm rate as a function of log 10 (L) when the 
target is described by a 2 x 2 template {N = K = 4). The threashold 7 is determined 
from 3.2 when the target template false alarm rate 1 satisfies P F (Target Template) = 

_ 1 Note that if the target template we^ a single pixel, the image false alarm rate would be approxi- 
mately 0.1 for all image sizes when Pf(Target Template) = Q.I/L-. 
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0.1/L 2 . The results shown in Figure 3.1 are obtained with 10,000 Monte-Carlo simu- 
lations. Note from Figure 3.1 that the image false alarm rate increases almost linearly 
with log 10 (L) and approaches 0.1 when L 2, i.e., when the image size is much larger 
than the template size. 

Figure 3.2 is similar to Figure 3.1 except that P F (Target Template) is ten times 
smaller, i.e., P F (Target Template) = 0.01/I 2 . The results shown in Figure 3.2 are 
obtained with 100,000 Monte-Carlo simulations. Note that a similar relation exists 
between the image false alarm rate and the target template false alarm rate. 

Figure 3.3 is also similar to Figure 3.1 except that the target is described by a 
3x3 template (N = K = 9). Note again that a similar relation exists between the 
image false alarm rate and the target template false alarm rate. 

The implication we obtain from the above three simple simulation examples is 
that when the image size is much larger than the target template size and a constant 
image false alarm rate P F (Image) (P F (Image) <C 1) is desired, the target template 
false alarm rate P F (Target Template) may simply be chosen as P F (Image) divided 
by the number of pixels in the image. The so obtained image false alarm rate will 
be approximately the same as or slightly smaller than the desired image false alarm 
rate. 
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Figure 3.1: Image false alarm rate as a function of log 10 (.L) when the target template 
false alarm rate is 0.1 [L 2 . The image size is Lx L and the target template size is 
2x2. The symbols denote the Monte-Carlo simulation results. The solid line 
denotes the straight line fitting, which has the form 0.0372 log 10 (X) + 0.0179, of the 




Figure 3.2: Image false alarm rate as a function of log 10 (Z,) when the target template 
false alarm rate is 0.01 /L 2 . The image size is L x L and the target template size is 
2x2. The symbols denote the Monte-Carlo simulation results. The solid line 
denotes the straight line fitting, which has the form 0.0058 log 10 (L) + 0.0009, of the 




Figure 3.3: Image false alarm rate as a function of log 10 (I) when the target template 
false alarm rate is 0.1/T 2 . The image size is L x L and the target template size is 
3x3. The “*” symbols denote the Monte-Carlo simulation results. The solid line 
denotes the straight line fitting, which has the form 0.0448 log 10 (L) - 0.0078, of the 
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4. An Efficient Algorithm for Two-Dimensional Frequency Estimation 


4.1 Introduction 

The two-dimensional (2-D) frequency/angle estimation algorithms, which have 
been studied recently, include the minimum variance method [1], the linear prediction 
methods [2, 3, 4, 5], the subspace rotation methods [6, 7, 8, 9, 10, 11, 12, 13, 14], 
the nonlinear least squares fitting method [15], and the maximum likelihood (ML) 
methods [16, 17]. These algorithms consider one of the following two data models. 
The first data model assumes that the frequencies of the 2-D complex sinusoids occur 
at arbitrary points of the 2-D frequency domain and the data model has the form 

= (4.1) 

£=i 

where {a^} denote the complex amplitudes and {ti^} and {/q:} denote the 2-D fre- 
quencies. For this data model, the number of unknowns is AK. The second data 
model assumes that the frequencies of the 2-D complex sinusoids occur at the inter- 
sections of some unknown grid lines in the 2-D frequency domain and the data model 
has the form _ _ 

= tt ( 4 -2) 

* =1 k=l 

For this data model, the number of unknowns is K+K +2/i K . Depending on the true 
distribution of the frequencies of the 2-D complex sinusoids, one model may result 
in fewer unknowns than the other, which, according to the parsimony principle [18], 
makes one model yield more accurate frequency estimates than the other. Consider an 
example where the frequency locations of the complex sinusoids are shown in Figure 
4.1. The numbers of unknowns in the first data model for the cases in Figures 4.1(a) 
- (g) are 12, 16, 20, 24, 28, 32, and 36, respectively. The numbers of unknowns in 
the second data model for these cases are the same and are 24. 


50 




(a) (b) <«0 < d > 




(e) (f) (8) 

• location of a sinusoid 


Figure 4.1: Locations of sinusoids in the two-dimensional frequency domain. 

This paper considers the second data model, which is also considered in [3, 5, 13, 
17]. This data model may be used when a vertical 2-D rectangular uniform linear 
array is used to estimate the incident angle of a signal arriving from a low angle 
relative to a smooth reflecting surface [19]. For this case, the signals arriving at 
the array consist of both the original incident signal and the signal reflected from 
the smooth surface. The same data model may also be used in synthetic aperture 
radar (SAR) or inverse SAR (ISAR) imaging [20, 21] to estimate the locations of the 
scattering centers of an object being imaged [5, 13]. These two cases will be further 
addressed in the paper. 

Among the previously studied methods, the subspace rotation methods, which 
include the state-space [6, 7, 9], ESPRIT [8, 10], and matrix pencil methods [11, 
13, 14], are known to be computationally efficient and have high resolution. Yet 
their estimation performances are usually not as good as that of the ML methods. 
The ML methods are asymptotically statistically efficient [17], but they require a 
computationally intensive multidimensional search over the parameter space. 
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This paper presents a computationally efficient eigenstructure-based 2D-MODE 
algorithm for estimating the frequencies of 2-D complex sinusoids. Unlike the sub- 
space rotation methods, the algorithm is statistically efficient under either the as- 
sumption that the number of temporal snapshots is large or the signal-to-noise ratio 
(SNR) is high. The statistical performance of the 2D-MODE algorithm is better than 
that of the computationally efficient subspace rotation methods, whereas the amount 
of computations required by the 2D-M0DE algorithm is usually no more than a few 

times of that needed by the subspace rotation methods. 

The remaining of this paper is organized as follows. In Section 2, we describe the 
problem of interest. In Section 3, we present the computationally efficient 2D-MODE 
estimator. In Section 4, we describe the properties and the applications of the 2D- 
MODE estimator. In Section 5, we establish the asymptotic statistical performance 
of the estimator under either the assumption that the number of temporal snapshots 
is large or the SNR is high. In Section 6, we provide several numerical examples 
comparing the performance of the 2D-M0DE algorithm with that of the subspace 
rotation algorithms. To make the paper self-contained, the subspace rotation algo- 
rithms we use are briefly described in the Appendix. Finally, Section 7 contains our 

conclusions. 


4.2 Problem Formulation 

Consider the following model of 2-D complex sinusoids in additive noise: 

„ K (i „) -tt ( 4 - 3 ) 


Vn 


where m = 1,2 m = 1,2 , ■ • • , M, and n = 1,2, •••,*■ We refer to M 
(M > I\) and A T ( M > A) as the numbers of spatial measurements , and to N as 
the number of temporal snapshots. The additive noise e m , w {t n ) is assumed to be a 
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complex Gaussian random process with zero-mean and 

£{e m ,m(*nJC, m^nj} =cr2 £n,,n 2 , (4-4) 

where (•)* denotes the complex conjugate and 6 ni ,„ 2 denotes the Kronecker delta. The 
em,m(f»), m = 1,2, • • • , M, rn — 1, 2, • • • , M, are assumed to be independent of each 
other and the complex sinusoids. The complex amplitudes a k,k(.^ n )'< k = L 2, ■ ■ • , A , 
k - 1,2, • • ■ ,7T, may be modeled as zero-mean complex Gaussian random processes 
that may or may not be correlated with each other. This model is referred to as the 
stochastic (or unconditional) signal model [22, 23]. For the stochastic signal model, 
the temporal signal snapshots at different sampling times t n are assumed to be 

independent of each other. The complex amplitudes may also be modeled as 

deterministic unknowns. This model is referred to as the deterministic (or conditional) 
signal model [22, 23]. 

Let 


y {tn) — [ yi,l{tn) ' ’ ' J/l,A/(^ n ) 


(4.5) 

e(fn) = [ ei ? i (tn) *•' e l,A ?(^n) 


(4.6) 

S (^n) = [ Q!i j i(f n ) <Xi t K(tn) 


(4.7) 

A = [ a(u?i) a(u>jc) 1? a ( u; fc) 

= [ e jwk ■ ■ ■ e jMuJk ] T , * = 1,2,-- 

■,K, 


(4.8) 


and 

B = [ b(nx) b(/tjr) ], b(pjf) = [ e 3 ^ e ]M ^ } T , I = 1 , 2 , ■ • • , A , 

(4.9) 

where (-) T denotes the transpose. Then ((4.3)) may be rewritten as 

y(f„) = (A®B)s(i n ) + e(<„). (4.10) 
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where ® denotes the Kronecker product. 

The problem of interest herein is to estimate u>i, u> 2 , U K an d ft 1 , ft 2 , ' ' 
from y (t n ), n = 1,2, Once the frequency estimates are obtained, the 

complex amplitudes of the 2-D complex sinusoids can also be estimated with a simple 
least-squares method [3]. 

4.3 The Computationally Efficient 2D-MODE Estimator 

The exact ML estimator of the 2-D frequencies requires a multidimensional search 
over the parameter space [17]. Since the search over the parameter space is compu- 
tationally prohibitive, computationally more efficient methods are attractive. We 
describe below a computationally efficient eigenstructure-based technique, which is 
an extension of a 1-D algorithm known as MODE [24] or, in a related form, WSF 
[25]. This estimator is herein referred to as the 2D-MODE. 

Let N 

R=i^ y (i„)y*(i„), (4.H) 

n=l 

where (-) H denotes the complex conjugate transpose and R is the estimate of the 
spatial covariance matrix R, 

R = E{y(t„)y H {t n )}- ( 4 - 12 ) 

It has been shown in [24, 25, 23, 26] that an asymptotically (for large N or high 
SNR) statistically efficient estimator of the 2-D frequencies u = [wi, wa, ■ ’ • ,“ K \ T 
and n = [fii, H 2 , • • ■ , ^] T can be o btained b y minimizing the following function: 

f(u>, fi) = tr [Pa 0B ( w . A*)E,A^A _1 Ef] , (4.13) 

where, for some matrix X, the symbol P£ stands for the orthogonal projector onto 
the null space of X, and where the columns in E s are the signal subspace eigenvectors 
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of R that correspond to the I\ largest eigenvalues of R, with A defined as 

K = min[iV, rank(S)]. (4.14) 

Here S is the signal covariance matrix, 

S = E{s{t n )s H (t n )}. (4-15) 

We assume that K is known. (If K is unknown, it can be estimated from the data as 
described, for example, in [27].) Note that if no components of the signal vector s(t) 
are fully correlated to one another and if N is large enough, then A = KK . Further, 

in ((4.13)) A is a diagonal matrix with diagonal elements Ai > A 2 > • • • > which 

are the I< largest eigenvalues of R, and 

A s = A — a 2 !, (4-16) 


with I denoting the identity matrix and 


MM - 1< i= %, 


MM - 1< 


It is worth noting that the evaluation of E*, A, and k 3 introduced above requires 
only the computation of the I< principal eigenpairs of R. Since usually I< < MM, 
the involved computational burden is of the order 0(M 2 M 2 ) and hence much reduced 
compared with what would be required for a full eigendecomposition. The orthogonal 
projector P^ 0B ma y be written as ( see f ° r the P ro P ertieS ° f the Kronecker 
product) 


Pa s b = I-(A®B)[(A ff ®B ff )(A®B)]“ 1 (A^®B ff ) 

= I — [a(a h a) -1 a h ] ® [B(B f B) _i B h ] 

= I-(I-Pa)®( i - p b) 

= I® P B + p a ® 1 _ p a ® Pb- (4-18) 
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Thus ((4.13)) may be rewritten as 


/ = tr [(i 0 P£) E, AjA-'Ef] +tr [(?i 0 i) E.A’A-'E"] -tr [(Pi ® P&) E, AJA 'E" 


(4.19) 

It is shown in Appendix A that the third term above is a higher-order term for 
both the case of large N and the case of high SNR. Thus minimizing ((4.19)) is 
asymptotically equivalent to minimizing 


/ = tr [(I ® P^) E,A?A-'E? ] + tr [(Pi 0 1) E s A^A -1 Ef] . (4.20) 


We show below how to avoid the search over the parameter space needed to 
minimize ((4.20)). The projector above may be reparameterized in terms of the 
coefficients c = [ Co _ Cl • ■ • c K V of a polynomial defined as 


y ft**’* = c„ n (z - e’"*) ; c, / 0. 

k= 0 k-\ 

Let C H be the following (M — K ) x M matrix 



0 


c 1 


CK 


Co 0 

Cl Co 


(4.21) 


(4.22) 


Then = Pc, where Pc denotes the orthogonal projector onto the range space of 
C (see [24, 29]). By defining d and D for B similarly as c and C for A, we obtain 
Pb = Pd- Thus by using the reparameterization above, ((4.20)) may be rewritten 
as 


/ = tr{[!0 (D(D h D)- 1 D b )] E,A;A-Ef)+tr{[(C(C''cr 1 C ,f ) ®I 


ELA^A _1 E k 


}• 


(4.23) 


Since E s is a consistent estimate of the true signal eigenvector matrix E s and 
the columns of E s belong to the range space of A ® B, it can be shown [24] that. 
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minimizing ((4.23)) is asymptotically equivalent to minimizing 

f = tr { [I ® (DW;'D')j E.AjA-'Ef } + tr { [(CWf'C*) ® i] E.A^A-'E" } , 

(4.24) 

where W x and W 2 are consistent estimates of C H C and D ff D, respectively. The 2D- 
MODE algorithm minimizes ((4.24)), in the following two steps. In the first step, let 
Wi = I and W 2 = I in ((4.24)) and minimize ((4.24)) to obtain consistent estimates 
c and d of c and d [29]. In the second step, let W x = C H C and W 2 = D H D, 
where C and D are formed from the c and d obtained in the first step, and minimize 
((4.24)) again. The c and d obtained in the second step are then used to compute 
the 2-D frequency estimates by rooting the corresponding polynomials of the form of 
((4.21)). 

Let j- 

Z„ ■ ■ ■ Z ijw 

Z = E,A2A-'Ef= ! ■■■ ! , (4.25) 

Z Ml • ■ ■ Z mm 

where Z,j are M x M matrices. Then 

tr(Zn) tr(Z x m) 

(CW^C^) : : > = nun ^{(CW^C*) Z w }, 

tr(Zjwi) tr(Z mm) 

(4.26) 

and 

d = mjn tr|(DW 2 -'D K ) j ^ mm tr {(dW X 'D") Z M } . (4.27) 

As shown in Appendix B, as either N or SNR increases, Z^ approaches a matrix Z ^ 
whose rank is no more than R. However, the rank of Z^ may be more than K. To 
improve the accuracy of the estimates of /i, we may replace Z ^ with its nearest rank- 
K approximation in the Frobenius norm metric [30]. Similarly, Z^j may be replaced 
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with its nearest rank- A" approximation in the Frobenius norm metric since the rank 
of the limiting matrix Z w is no more than K (as also shown in Appendix B). For 
detailed implementation procedures of minimizing the quadratic functions ((4.26)) 
and ((4.27)), see [29, 31]. 

The polynomial coefficients in ((4.26)) and ((4.27)) can be constrained by making 
use of the so-called conjugate symmetry conditions [29, 31] 

c k = c* K _ k , k = 0,1, (4-28) 

and 

<% = <%-!. * = ( 4 - 29) 
These constraints, considered in detail in the cited works, leads to the most parsi- 
monious parameterization of the estimation problem under discussion, and hence can 
yield enhanced estimation performance [18]. In the simulation examples given in Sec- 
tion 6, we present results obtained by imposing the conjugate symmetry constraints 
on the polynomial coefficients in ((4.26)) and ((4.27)). 

4.4 Applications and Properties of the 2D-MODE Estimator 

We consider below the use and the properties of the 2D-MODE estimator in the 
case where the number of temporal snapshots N is large and the case where the SNR 
is high. Since in practical applications, particularly those involving sensor arrays, it 
is less likely to have large M or M, we shall not consider such a case herein. 

4.4.1 The Case of Large N 

The case where the number of temporal snapshots N is “large” may occur when 
estimating 2-D incident angles with a 2-D MxM rectangular uniform linear array. For 
this application, the {u k } and { n } are the phase factors from which we can calculate 
the 2-D incident angles (i.e., the azimuth and elevation angles). Equation ((4.3)) 
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models certain special cases where the incident angles are related to each other. For 
example, such a special case occurs when a vertical 2-D rectangular uniform linear 
array is used to estimate the incident angles of a signal arriving from a low angle 
relative to a smooth reflecting surface such as the calm sea [19). For this case, the 
signals arriving at the array consist of both the original incident signal and the signal 
reflected from the surface, i.e., the specular path, as shown in Figure 4.2. The incident 
angle and its reflected angle are related to each other and can be shown to fit into 
the data model in ((4.3)). For this application, N denotes the number of independent 
snapshots taken at the output of the array. Furthermore, 

Uk = 9k sin 4>k, (4.30) 


and 


n = 


2tt8i 

^0 


cos 6 k , 


(4.31) 


where \ 0 denotes the wavelength of the incident signals and 8 k and ds denote the 
elevation and azimuth angles, respectively. Since di = h and <h = 180"-»,. we have 
K = 1, H = 2, and pi = -H for this example. 

We can also use the knowledge that H = -/ij to obtain improved estimates when 
minimizing ((4.27)). With pi = -/l 2 , we can constrain the polynomial coefficients d, 
in ((4.27)) so that * = d, = 1 and d, is real. Note that these conditions are more 
strict than the conjugate symmetry condition in ((4.29)). Under these conditions, the 
minimization of ((4.27)) may be carried out as follows. Let ZpZp denote the rank-2 
approximation of Z„, where Z M is an M x 2 matrix. Let i ! P and zP be the first and 
second columns of Z/z, respectively. Then for k = 1,2, 



■ .IX . ~n 

7 - 4- z — 
Z 3 ,k ^ \,k 

-H 1 

<fc 






zl* 

*4,1 + Z 2,k 

5^- 

3,k 


1 

di 1 

=H-r' 

do 


~M,k ^ M-2,k 

Z ~M - 1 ,fc - 

1 




(4.32) 
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where 5^, i = 1, 2, • • ■ , M, denotes the zth element of . Let 


Re{W^ /2 Hf} 

Re{W^ /2 Hf} 

Im{W^ /2 Hf} 

Im{W^ /2 Hf} 


-ft 

gi 



(4.33) 


where j = 1, 2, denotes the jth column of G^. Then minimizing ((4.27)) becomes 


equivalent to minimizing 
the solution is 




1 

d\ 


where II ■ || denotes the Euclidean norm, and 


d , = - [(gf)"^]" 


(4.34) 


In the case where the number of temporal snapshots N is large, similar to the 
results obtained in [24, 23] for the 1-D angle/frequency estimation case, the compu- 
tationally efficient 2D-M0DE estimator can be shown to be statistically equivalent to 
the stochastic maximum likelihood (ML) estimator (i.e, the ML estimator obtained 
under the stochastic signal model). Since the stochastic ML estimator asymptotically 
achieves the stochastic Cramer-Rao bound (CRB), it follows that the 2D-M0DE es- 
timator is an asymptotically (for large N) statistically efficient estimator. 

We remark that the assumption that the temporal signal snapshots a k g(t n ) are 
independent of each other is not a necessary condition for the application of the 
2D-M0DE algorithm. In particular, the asymptotic accuracy of the 2D-M0DE algo- 
rithm is the same for both temporally correlated and temporally uncorrelated signal 
snapshots (see Remark 3 on p. 1787 of [22]). For the temporally correlated sig- 
nals, however, more accurate estimators may be devised if the temporal correlation 
function is known or if it is known to have some parameterized functional form. 
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4.4.2 The Case of High SNR 

The case where the SNR is “high” and N is small may occur in synthetic aperture 
radar (SAR) or inverse SAR (ISAR) imaging [20, 21, 13, 5]. In this application, the 
radar usually transmits linear frequency modulated (chirp) pulses. Upon receiving 
each pulse returned by an object being imaged, the radar mixes the pulse with a 
reference chirp signal and low-pass filter the mixed signal. As a result, the scattering 
centers of the object at different ranges correspond to different frequencies of the 
output of this operation. Since either the radar or the object is moving or rotating, 
the pulses received at different angles between the radar and the object are used to 
form a synthetic aperture. The scattering centers of the object at the same range but 
different cross-ranges correspond to different (Doppler) frequencies over the synthetic 
aperture. Thus for this application, the {u k } and {/z*} describe the locations, i.e., 
the ranges and cross-ranges, of the scattering centers of the object being imaged. 

Similar to the results obtained in [26] for the 1-D angle/frequency estimation case, 
the computationally efficient 2D-MODE estimator can be shown to be asymptotically 
(for high SNR) statistically equivalent to the deterministic ML estimator. It has 
been shown in [26] that for high SNR, the deterministic ML estimator achieves the 
deterministic CRB. Thus the 2D-MODE estimator is an asymptotically (for high 
SNR) statistically efficient estimator in this case also. 

We remark that in the case of high SNR, the N may take any value greater than or 
equal to 1. The case of N = 1 has a particular relevance for the SAR/ISAR application 
mentioned previously [20, 21]. In such a case, we have E 5 = y(^i)/||y(U)|| &nd A 5 
and A are scalars. Thus for N = 1, minimizing ((4.13)) is equivalent to minimizing 

/ = y H (*i)PA®By(*0- ( 4 - 35 ) 

We note that this function is exactly the one to be minimized by the deterministic 
ML estimator given in [17]. Since the 2D-MODE algorithm determines an asymptotic 
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realization of the minimizer of ((4.13)), it readily follows that the 2D-MODE and the 
deterministic ML estimators are asymptotically (for high SNR) equivalent. (We note 
that showing this equivalence in the case of N > 1 is a more complicated operation 
[26].) 

4.4-3 Parameter Identifiability Conditions 

We note that our previous results are valid if M and M are constrained by the 
parameter identifiability requirements [23, 32, 33]. Let Tj denote a real- valued vector 
containing all unknowns of the data model in ((4.3)). Then our problem is parameter 


identifiable if 


R — Ro ^ k T] 77 0 . 

(4.36) 

As shown in Appendix C, the inequalities 


M > 2 K - =, 
M 

(4.37) 

and 

M > 2K -f-, 

(4.38) 

are sufficient conditions for parameter identifiability. 

By rearranging the elements of y (f n )? s(£ n ), and e(£ n ), we can 
y(t n ) == (B ® A)s(£ n ) + e(f n ). We can then show similarly that 

rewrite ((4.10)) as 

M> 27- A, 
M 

(4.39) 

and 

M > 2 K - ==, 

(4.40) 


are also sufficient conditions for parameter identifiability. 
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4.5 Statistical Performance Analysis 

In this section, we establish the asymptotic statistical performance of the 2D- 
MODE estimator for both the case of large TV and the case of high SNR. 

4.5.1 The Case of Large N 

As argued in Appendix D, the asymptotic (for large TV) statistical distribution 
of Cj is Gaussian with mean w and covariance matrix equal to the corresponding 
stochastic Cramer-Rao bound (CRB), CRB^,. The ijth element of (CRB^) 1 is 
given by 

[(CRB^,)- 1 ] .. = |^Re (tr { [(Af Pi A,) ® (b"b)] S (a" ® B") R- (A ® B) S}) , 

(4.41) 

where A,- = 5A Similarly, the asymptotic (for large N) statistical distribution 
of jx is shown to be Gaussian with mean fx and covariance matrix equal to the corre- 
sponding stochastic CRB, CRB^. The ijth element of (CRB^) 1 is (see Appendix 

D) 

[(CRBy-’].. = |^Re(tr{[(A»A) ® (BfPjjB,)] S (a» ®B h ) R- (A® B)s}) , 
U C (4.42) 

where B; = <9B/d/x t . 

4.5.2 The Case of High SNR 

It is shown in Appendix D that the asymptotic (for high SNR) statistical distribu- 
tion of Cj is Gaussian with mean uj and covariance matrix equal to the deterministic 
CRB given by 

(CRBi,)" 1 ] = — -^-Re {tr ([(AgPiA fcl ) <2> (B h b)] s)} , (4.43) 

. ' ' J i j & 
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where 


(4.44) 


S = ^ E s(<n)s W (<n) 

Tl=l 

It is also shown in Appendix D that the asymptotic (for high SNR) statistical distribu- 
tion of is Gaussian with mean p and covariance matrix equal to the corresponding 
deterministic CRB given by 

'(CRB-,)' 1 ] = ^Re(tr{[(A»A) ® (BfPiBi)] S}) . (4.45) 

4.6 Numerical Examples 

In this section, we illustrate the performance achievable by the 2D-M0DE algo- 
rithm and compare it with the performance of the two-dimensional subspace rotation 
methods (2D-SRMs). (The 2D-SRMs, both with and without spatial smoothing, are 
briefly described in Appendix E.) The empirical performance of the algorithms is ob- 
tained from 100 independent trials and it is compared with the theoretical statistical 
performance given by the CRBs. 

In the first three examples, we consider the case of large N . In particular, we 
consider the scenario shown in Figure 4.2 where a vertical rectangular uniform linear 
array with M = 8 and M = 10 is used to estimate the angles of arrival of an incident 
signal arriving from (<M) and of its reflected signal arriving from (<f>, 180° - 9). The 
signal-to-noise ratio (SNR) of the direct signal used in the examples is assumed to be 
-7 dB. The SNR of the reflected signal is assumed to be 3dB less than that of the 
direct signal. The spacings and S 2 between two adjacent sensors in the array are 
assumed to be a half wavelength. The CRBs of the estimates of <f> and 9 are readily 
obtained from the CRBs of the estimates of u k and K given in Section 5 and the 
Equations ((4.30)) and ((4.31)) relating <f> and 9 to u k and /xp 

Figure 4.3 shows the root-mean-squared errors (RMSEs) of the angle estimates 
as a function of the elevation angle 9 when (f> = 45° and N = 60. The incident and 
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reflected signals are assumed to be uncorrelated with each other, which may happen 
when the rectangular array is well above the reflecting surface. For this example, 
the 2D-SRM is used without spatial smoothing. We note that the performance of 
the 2D-M0DE is always better than that of the 2D-SRM. Figure 4.3 also shows the 
asymptotic (large JV) statistical performance, corresponding to the CR-bound (CRB). 
We note that the 2D-M0DE estimator performance is very close to its asymptotic 
statistical performance, which is also the best achievable performance in the class of 
(asymptotically) unbiased estimators. 

Figure 4.4 shows the performance of the 2-D estimators as a function of the cor- 
relation coefficient between the incident and reflected signals when <f> = 45°, 0 = 85 , 
and N = 500. For this example, the 2D-SRMs are used both with and without 
spatial smoothing. We note that the 2D-MODE algorithm has the most significant 
advantage over the 2D-SRM without spatial smoothing when the signals are highly 
correlated or coherent. When the incident signals are completely correlated with 
each other or coherent, the 2D-SRM without spatial smoothing fails. For this case, 
the 2D-SRM must be used with spatial smoothing. For the 2D-SRM with spatial 
smoothing, we set both dimensions of the subarrays to L = L = 5. Figure 4.4 shows 
that the 2D-SRM with spatial smoothing is also outperformed by the 2D-MODE. In 
our implementations, the amount of computations needed by the 2D-MODE is about 
7.1 times of that needed by the 2D-SRM with or without spatial smoothing. 

Figure 4.5 shows the RMSEs of the angle estimates as a function of N when 
<j) = 45°, 6 = 85°, and the incident and reflected signals are 99% correlated. For 
the 2D-SRM with spatial smoothing, both dimensions of the subarrays are set to 
L = L = 5. We note again that the 2D-MODE algorithm performs better than the 
2D-SRM both with and without spatial smoothing. We also note that the larger the 
N , the closer the 2D-MODE estimator performance to its CRB. Note from Figures 4.3 
and 4.5 that for the 2D-MODE estimator to achieve its CRB, a larger N is needed for 
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h.ghly correlated signals than for nncorrelated signals. For N = 10,000, the amount 

o co Ut by the 2d mode . about 4 o timK of that nMdsd 

2D-SRM with or without spatial smoothing. 

la the next three examples, we consider the performance of the 2D- MODE fre 
quency estimator for the case «T * - which is the radar imaging case. We shall 

r/ h r r are tw * • - - - - two 

assume that M = 8 and^ - ^ ^ "" ^ ** *° 

, , 10 ' The P erf °™ance of the 2 D-MODE estimator is 

compare wit that of the 2D-SRM with spatial smoothing (with L = I = 5 ). Note 

tha as for the case of coherent incident signals and large the 2D-SRM must be 

used with spatia, smoothing for this = , case, no matter what the SNR is 

Figure 4.6 shows the RMSEs of the frequency estimates as a function of SNR 

that t Tf = X 0 ' 28 ’ X °' 33) an<1 M = <2T X °' 30 ' ^ note 

a «h performance of the 2D-MODE algorithm is a, ways better than that of the 

, “ SPaHal Sm ° 0thing ’ eSPeCid,y W, “ SNR * -und 0 dB. Figure 4 6 

a^ shows the asymptotic (high SNR) statistical performance, corresponding to the 

it ‘ °“ (CRB) - ^ a ° te tha ‘ the 2DM0DE -‘in-tor performance is very dose to 
s asymptorfc statistical performance, which is also the best achievable performance 

or “ eStimatOTS ’ “• SNR * -ater than 

Figure 4.7 shows the RMSEs of the frequency estimates as a function of the fre 
quency separation A w /(2x) when (w,,^) = (2r x 0 .2S,2g x 0.28 + Aw) (p, _ 

(2ir x 0.30, 2x x 0.35), and SNR = 10 dB Figure 4 8 sh ,i , ’ 

as a fund- rr ae. figure 4.8 shows the frequency estimates 

function of frequency separation Ap/(2x) when („, , w 2 ) = (2r x 0 . 2S x 0 33) 

li.:i ; ?: »— ** - — »» ».L ^^2 

t e -D MODE algorithm again is always better than that of the 2D-SRM 
SPatial Sm °° thinS ’ ^ frequency separation, either ^ 
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or Ap/(2tt), is small. We also note that as the frequency separation increases, the 
2D-M0DE estimator performance becomes closer to its asymptotic statistical perfor- 

mance. 

Because the eigendecomposition of R is no longer needed for N = 1, the amount 
of computations needed by the 2D-MODE for N = 1 is reduced to 0.35% of the 
amount required by the example in Figure 4.4. The amount of computations needed 
by the 2D-SRM with spatial smoothing for N = 1 is also reduced to 0.80% of the 
amount required by the example in Figure 4.4. For N = 1, L = L = 5, M = 8, and 
M = 10, the amount of computations needed by the 2D-MODE is about 3.1 times of 
that needed by the 2D-SRM with spatial smoothing. 

4.7 Conclusions 

We have presented a computationally efficient eigenstructure-based 2D-MODE 
algorithm for two-dimensional frequency estimation. We have shown that this esti- 
mator is asymptotically statistically efficient under the assumption that either the 
number of temporal snapshots is large or the signal-to-noise ratio is high. Numerical 
examples showing the comparative performances of this algorithm and of the compu- 
tationally efficient subspace rotation algorithms have also been given in this paper. 
We have shown that the performance of the 2D-MODE algorithm is better than that 
of the subspace rotation methods, whereas the amount of computations required by 
the 2D-MODE algorithm is usually no more than a few times of that needed by the 
subspace rotation methods. 

Appendix A - The Simplification of the 2D-MODE Cost Function 

We show below that the third term of ((4.19)) is a higher-order term and may be 
neglected for large A or high SNR. 
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The Case of Large N 


Consider first the case of N > 1. Let /i, / 2 , and [:> denote the first, second, and 
third terms, respectively, of the right side of ((4.19)). Let (■)' denote the gradient of 
(•) with respect to fJ. Since A minimizes ((4.19)), we have for N > 1. 

A - M = -(/TV' = + /T'(/; + &)■ < 4 -«) 

The zth element of /{ is 

1/iK = tr{(l® [Pfe]') E.A^A-'E"} , (4.47) 

where [25] 

[p^]' = -B tH [B ,? l'P^-(B w |B H ];Py' f , (4.48) 

with 

B r =(B H B) _1 B. (4-49) 

Thus [fi\i may be rewritten as 

= -2 Re (tr { [i ® (B tH [B H ][PB)] E s A*A _1 Ef }) . (4.50) 

Since the columns of E, span the same signal subspace as the columns of A ® B and 
PgB == 0, we have (i ® Pfj) E s = 0. Thus 

M = -2 Re (tr { [i <g> (B^B*]^)] (E. - E.) A^A _1 Ef}) . (4.51) 

Since (E s - E s ) = 0(1/ y/N) [34, 35], we have [/i]- = 0(l/\/A"). The zth element of 
f’ z may be written as 

[/ 3 ]' = -2 Re (tr { [Pi <g> (B tH [B H ]'P£)] (E 3 - E s ) A* A' 1 (E s - E s ) j) , 

V 1 (4.52) 
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where we have used the fact that (P^ ® i) E, = 0. Thus [/ 3 ]' = 0(1 jN) is negligible 
for TV > 1 as compared with [/i]'-. It can be shown similarly that [/ 3 ]"- is also 
negligible for N 1 as compared with [/i];p 
Thus for N 1, 

a- a *=- c/rrvi, ( 4 - 53 ) 

which means that the third term of ((4.19)) is a higher-order term when estimating 
fj, and may be neglected asymptotically (for N » 1). We can show similarly that 
this result is also true when estimating w. 

The Case of High SNR 

Consider next the case of high SNR, i.e., a < 1, where a is the standard deviation 
of the additive noise. For a <C 1, Equation ((4.46)) also holds. Since (i® Pb) — 
(I ® Pfe) (E. - E.) = 0(a) and (P* ® i) E, = (Pi ® i) (E, - E.) = 0(a), we 
obtain [h]' = 0(a) from Equation ((4.51)) and [/ 3 ]' = 0(a 2 ) from Equation ((4.52)). 
Thus [/ 3 ]' is negligible for cr < 1 as compared with [/i]'. It can be shown similarly 
that [/ 3 ]"- is also negligible for a < 1 as compared with [/i]-'-. Then the third 
term of ((4.19)) is a higher-order term when estimating n and may be neglected 
asymptotically (for a < 1). It can be shown similarly that this result is also true 

when estimating w. 

Appendix B — The Ranks of Zu; and Z^. 

To determine the ranks of Zu; and Z^t, let us consider 

Zn * ■ * Z xM 

Z = E s AjA _1 E^ = i : 

Z.\/i • ■ ■ Timm 


(4.54) 
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where Z,y are M x M submatrices of Z. Note first that E s may be written as 

E s = (A<g>B)Q, (4-55) 


where Q is some (KI<) x K matrix of full rank K = min[iV,rank(S)] with probability 
1. Thus Z may be rewritten as 

Z = {A®B)QA 2 s A- 1 Q h {A h ®B h ). (4.56) 


Consider first the rank of Z^. Note that 

Ban 

I 

(A0B) = 


Baift- 


(4.57) 


Bc&a/i ■ ■ ' B cimk 

where a,, denotes the ijth element of A. Thus Z„ has the form B(- • -)B" and 


Z M =i;Z„„=B(..-)B H . (4.58) 

m = 1 

Then the rank of Z^ is no more than the rank of B. 

Consider next the rank of Zu;. Using (A ® B) = (A ®I)(I® B), we get 

Z = (A<E>I)r(A F <g>I), (4-59) 


where 


r = (1 0 B)QA^A _ 1 Q ff (1 0 B H ) = 


where T ,-j are M x M sub mat rices of T. Thus 
tr(Z,j) = tr { (ai r) 0 i) T [(a< r) ) H 0 1] } = ^ 


r n 

• T lK 

Tki ■■ 

' . 3 

• T K k 

tr(r u ) • 

•• tr(TiA-) 


tr(IVi) tr(IVA-) 


(4.60) 



(4-61) 


where a-* 4 denotes the zth row vector of A. Hence Z u ; has the form Z u — A(- ■ -)A H . 
Then the rank of Z^ is no more than the rank of A. 

Appendix C - A Sufficient Condition for Parameter Identifiability 

We prove below the sufficiency of the parameter identifiability requirements ((4.37)) 
and ((4.38)). Our proof follows closely the one for 1-D angle/frequency estimation in 

[32, 33]. 

As in [32, 33], to prove the sufficiency of ((4.37)) and ((4.38)), we need to show 
that ((4.37)) and ((4.38)) are sufficient conditions for A = A 0 , B = B 0 , and Q = Qo 
to be the unique solution of 

(A <g> B)Q = (A 0 ® B 0 )Qo, ( 4 - 62 ) 

where Q is some (I<K) x K matrix of full rank K. The Equation ((4.62)) may be 
written as 

(A <g> I) (I ® B)Q = (A 0 ® I)(I ® B 0 )Qo- ( 4 - 63 ) 

Let 

(I®B)Q = 0, and (I®B 0 )Qo = Oo. ( 4 - 64 ) 

The O is of dimension (KM) x K and is of full rank K since M > K and B is a 
Vandermonde matrix. Then ((4.63)) may be written as 

(A® 1)0 = (A 0 ®I)fio- ( 4 - 65 ) 

Consider first a solution A and O of ((4.65)), which is such that the columns *i, 
i 2 , • - i c in A and A 0 coincide, where 0 < c < K. Since c < K, the solution A = A 0 
and O = O 0 is excluded. Let A 0 be the submatrix of A 0 without the columns i u i 2 , 
• - -, i e . Then (Ao®I) is of dimension (MM) x (KM-cM). Let Q 0 be the submatrix 
of flo without the rows that correspond to the columns of (A <g> I) which are missing 
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in (A 0 ® I). Let ft be the same as fi except that the rows of fio that are missing in 
tl 0 are subtracted from the corresponding rows of ft. Thus ((4.65)) can be written 


as 




(4.66) 


Let 

71 = dim | A ^ (A® I) ( Ao ® I) )} > (4-67) 

where dim{A'(-)} denotes the dimension of the null space of a matrix. Since A is a 
Vandermonde matrix, the vectors in A for any M different values of u> k are linearly 
independent. Thus 

7l = 2 KM -cM - min [MM, (2 KM - cM )] . (4.68) 



(4.69) 


where dim{77.(-)} denotes the dimension of the range space of a matrix. Our goal is 
to prove that ((4.66)) has no solution under ((4.37)), or equivalently 


7i < 72- 


(4.70) 


If MM > (2/1 M - c¥), then 71 = 0 and ((4.70)) follows since 72 > L If MM < 
(2 KM - c¥) , then 71 = 2KM-cM-MM and from ((4.64)), we have 72 > K-cM. 
From ((4.37)), we have 2 I<M - cM - MM < K - cM. Thus under ((4.37)), the 
unique solution to ((4.65)) is A = Ao and O = fio- 
Consider next the solution of 

(I® B)Q = (I® B 0 )Qo- ( 4 - 71 ) 

We can show similarly that under ((4.38)), the unique solution to ((4.71)) is B = B 0 
and Q 7 Qo- 
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Thus the proof is finished. 


Appendix D - Derivation of the Cramer-Rao Bounds 

We present the derivations of the Cramer-Rao Bounds (CRBs) for both stochastic 
and deterministic signal models. 

Stochastic CRB 

It follows from [22, 23, 24] that the asymptotic (for large N) statistical distribution 
of the parameter estimates [ jt T ] T obtained with the 2D-M0DE algorithm is 
Gaussian with mean [ u T fi T \ T and covariance matrix equal to the stochastic CRB, 

CRB S . Let 

A = A®B. ( 4 - 72 ) 

Using Equation (4.63) in [23], we obtain the ij th element of (CRB 5 ) 1 as 

[(CRB 5 ) -1 ] .. = [ tr (^fPiASA H R - Ms)] , (4.73) 

with Ai denoting the derivative of A with respect to the ith element of [ u T [i ] 

We first show that the CRBs for Cj and /t are decoupled. Let Ak denote the 
derivative of A with respect to the kth element of u>. Then 

A k = A* ® B, (4-74) 

where A k = dA/du k . Let Ai denote the derivative of A with respect to the kth. 
element of fi. Then 

A % = AOB^, (4.75) 

where B^ = dB/d^. Using ((4.18)), we obtain 

AZ Pi An = (a^ ® B h ) (I ® P£ + Pa © 1 - p A ® Pb) ( a ® B *) 

= (A? A) ® (B»P^Bp) + (AfPiA) ® (B»Bj) - (A?PiA) ® (B*PiB r ) 
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(4.76) 


= 0 , 

which shows that the CRB S is block-diagonal. 

Next note that 

/'PjA, = (Ag®B H ) (I®Pb + Pa® i - p a® p b)( a *. ® B ) 

= (AgPiA*.) 0 (b"b) . (4-77) 

Hence the (ki, A? 2 )th element of (CRBy) 1 is 

[(CRBy)' 1 ]^^ = ^Re{tr ([(A^PiA fcl ) <8> (b h b)] S^R“Ms)}. (4.78) 

Thus using ((4.72)), we complete the proof of ((4.41)). The proof of ((4.42)) is similar. 

Deterministic CRB 

It follows from [26] that for high SNR, the statistical distribution of the esti- 
mate [ Cj t } T obtained with the 2D-M0DE algorithm is Gaussian with mean 
[ <jj T f. i T ] T and covariance matrix equal to the deterministic CRB, CRB . Using 
Equation (4.68) in [23], we obtain the ij th element of (CRB^) 1 as 

(CRB' , )-'] j . = Ui;Re [s"(t„)>lf p iA-s(«»)]} • (4.79) 

Using ((4.76)), we can show the CRB d is also block-diagonal. Thus the CRBs for u 
and ii are decoupled for the deterministic signal model as well. 

Using ((4.77)), we get the (k u k 2 ) th element of (CRB^) -1 as 

[(CRB4,)-'] Mi = 4|ERe[s' , (l„)((Ag p iA il )®(B ff B))s(i„)]J 

= ^R»{tr([(A5PiA h )®(B'B)]S)}. (4.80) 

Thus ((4.43)) is proven. The proof of ((4.45)) is similar. 


75 



Appendix E - The 2D-SRM Algorithms 

We briefly describe below the two-dimensional subspace rotation methods (2D- 
SRMs) we used in the numerical examples. The methods are similar to the matrix 
pencil method detailed in [13], which was derived for N — 1. 

The 2D-SRM without Spatial Smoothing 

For uncorrelated or slightly correlated signals and large N , the following method, 
which we refer to as the 2D-SRM without Spatial Smoothing , may be used to estimate 

the 2-D frequencies. 

Let Y(t n ) be an M x M matrix whose (m,m)th element is To estimate 

uj, let 

Y w = y w J (Y w )* ] » ( 4 - 81 ) 

where J denotes the exchange matrix (with ones on the antidiagonal and zeros else- 
where) and 

Y w = Y(<i) Y(t 2 ) Y{t N ) ■ ( 4 - 82 ) 

Let the columns of the M x I< matrix U s be the left singular vectors of Y w that 
correspond to the largest singular values of Y w . Let U s i and U 3 2 denote the (M — 
x K submatrices of U s consisting of the first and last (A/ 1) rows of U 3 . Then 

the eigenvalues of {V%U * are the estimates of e JWl , • • ■, e^ K . From 

these eigenvalues, we obtain the estimate of un 

The steps of estimating n are similar to those of estimating w. The mam difference 
is that the Y(f„), n = 1 , 2 ,* • • 5 iV, in ((4.82)) above are replaced by Y T {t n ). The 
parameters M and K are also changed to M and K, respectively. 
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The 2D-SRM with Spatial Smoothing 


For highly correlated signals or small N, the following method, which we refer to 
as the 2D-SRM with Spatial Smoothing , may be used to estimate the 2-D frequencies. 

To estimate w, let 



where 


L ss 


Y„(U) Y ss {ii) 



(4.84) 


with 


Y„(in) 


Y„l(t„) Y ss 2 {tn) 



(4.85) 


and Y 3S i{t n ), i = 1, 2, • • ■ , M — L + 1, denoting the submatrix of Y(t n ) that consists 
of the ith to (* + L - l)th rows of Y{t n ). The remaining steps of the 2D-SRM with 
Spatial Smoothing are similar to those of the 2D-SRM without Spatial Smoothing. 

The steps of estimating /r are again similar to those of estimating u>. The mam 
difference is that the Y„,((„), i = 1, 2, ■ ■ ■ , M - L + 1, in ((4.85)) above are obtained 
as submatrices of Y T (t „). The parameters M, K, and L are also changed to M, K, 
and L, respectively. 
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(a) 



(b) 

Figure 4.3: Root-mean-squared errors of angle estimates as a function of the elevation 
angle 6 when the direct and reflected signals arrive from (45°, 0) and (45°, 180° — 0), 
respectively. The SNRs for the direct and incident signals are —7 dB and —10 dB, 
respectively. The signals are assumed to be uncorrelated with each other. Further, 
M — 8, ~M — 10, and N = 60. (a) For the estimates of the azimuth angle cj>. (b) For 
the estimates of the elevation angle 6. 
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Figure 4.4: Root-mean-squared errors of angle estimates as a function of the correla- 
tion coefficient between the direct and reflected signals when the direct and reflected 
signals arrive from (45°, 85°) and (45°, 95°), respectively. The SNRs for the direct 
and incident signals are — 7 dB and —10 dB, respectively. Further, M = 8, M = 10, 
L = L — 5 (for 2D-SRM with spatial smoothing), and N = 500. (a) For the estimates 
of the azimuth angle 4>. (b) For the estimates of the elevation angle 6. 
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Figure 4.5: Root-mean-squared errors of angle estimates as a function of the number 
of temporal snapshots N when the direct and reflected signals arrive from (45°, 85°) 
and (45°, 95°), respectively. The SNRs for the direct and incident signals are -7 dB 
and — 10 dB, respectively. The correlation coefficient between the direct and reflected 
signals is 0.99. Further, M = 8, M = 10, and L = L = 5 (for 2D-SRM with spatial 
smoothing), (a) For the estimates of the azimuth angle 4 > . (b) For the estimates of 
the elevation angle 6. 
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Figure 4.7: Root-mean-squared errors of frequency estimates as a function of 

frequency separation Au/2ir when (wi,^) = (27 r x 0.28, 27r x 0.28 -+^Ao;), 
(^ 1 ,^ 2 ) — (27 r x 0.30, 2tt x 0.35), SNR = 10 dB, M = 8, M = 10, L = L = 5 
(for 2D-SRM with spatial smoothing), and N = 1. (a) For the estimates of 
(b) For the estimates of ni/(2n). 
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Frequency Separation 
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Figure 4.8: Root-mean-squared errors of frequency estimates as a function 

of frequency separation A^j/ 27 r when (uq,u> 2 ) = x 0.28,2? r x ^0 . 33), 

(Huli 3 ) = (2? r x 0.30, 2 tt x 0.30 + A/z), SNR = 10 dB, M = 8, M = 10, L = L = 5 
(for 2D-SRM with spatial smoothing), and N = 1. (a) For the estimates of w 1 /(2tt). 
(b) For the estimates of ^i/(27r). 





5. High Resolution Range Signature Estimation and Synthetic 

Aperture Radar Imaging 


5.1 Introduction 

This paper proposes a robust parametric data model for estimating high resolution 
range signatures of radar targets and for forming high resolution synthetic aperture 
radar (SAR) images. For the range signature estimation, estimating the radar cross 
section (RCS) of the scattering center of a radar target at a certain range is modeled 
as estimating the amplitude and phase of a complex sinusoid with known frequency 
in unknown colored Gaussian noise. For the SAR imaging application, estimating the 
complex intensity of a pixel in an SAR image is modeled as estimating the amplitude 
and phase of a two-dimensional complex sinusoid with known frequency in unknown 
colored Gaussian noise. This new modeling approach models the thermal noise and 
the interferences from other scattering centers of a radar target or other pixels in a 
SAR image as unknown colored noise. 

This paper also presents a parameter estimation algorithm for the data model. 
The algorithm is referred to as the APES (Amplitude and Phase Estimation of a 
Sinusoid in unknown colored noise) algorithm. We shall describe how the APES 
algorithm can be used to estimate the range signatures and to form SAR images. 
The APES algorithm avoids the search over the parameter space and requires only 

simple matrix multiplications and matrix inverses. 

We will show that our modeling and estimation approach yields better resolution 
and lower sidelobes than the conventional nonparametric FFT (fast Fourier trans- 
form) method. We will also show that our approach is more robust than modeling 
the radar data as a certain number of complex sinusoids in noise and estimating the 
frequencies, amplitudes, and phases of the sinusoids with one of the best sinusoidal 
parameter estimation methods. We will present both numerical and experimental 
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examples comparing the performance of our approach with the two afore-mentioned 
approaches. 

In Section 2, we consider high resolution range signature estimation. We shall 
formulate the problem of interest, present the APES algorithm, and demonstrate 
the performance of our approach with both simulated and experimental data. (In 
Appendix A, we derive a computationally and asymptotically statistically efficient 
large sample maximum likelihood estimator for estimating the complex gain of a signal 
with known waveform and known steering vector in unknown colored Gaussian noise, 
which the APES estimator mimics.) In Section 3, we consider high resolution SAR 
imaging with the APES method. We shall formulate the problem of interest, describe 
how to use the APES method for SAR imaging, and demonstrate the performance 
of our approach with XPATCH (a computational electromagnetics software) data. 
Finally, Section 4 contains our conclusions. 

5.2 High Resolution Range Signature Estimation 

In this section, we propose a parametric data model and present an estimation 
algorithm for estimating high resolution range signatures of a radar target. High reso- 
lution range signatures are useful for many applications including the non-cooperative 
target identification (NCTI). 

5.2.1 Problem Formulation 

The range resolution of a radar is determined by the radar bandwidth. To achieve 
high resolution in range, the radar must transmit wideband pulses, which are often 
linear frequency modulated (chirp) pulses. Upon receiving each pulse returned by a 
radar target, the radar demodulates the pulse by mixing the pulse with a reference 
chirp signal and low-pass filter the mixed signal. As a result, the scattering centers 
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of a radar target at different ranges correspond to different frequencies of the output 
of this demodulation operation. 

Let z denote an Mx 1 vector containing the samples of the output of this operation. 
We assume that the radar bandwidth is moderate so that the radar cross sections 
(RCSs) of the scattering centers of the radar target do not change with frequency. 
Then z may be written as 

z = + n o, (5-1) 

k=i 

where 7 *, k = 1 , 2, • • • , I < , denotes the RCS of the kth scattering center of the radar 
target, no denotes the additive noise vector, and 


h 


1 e~ jtk 





(5.2) 


with (•) r denoting the transpose and t k denoting the time delay proportional to the 
range of the kth scattering center. 

A simple nonparametric method of estimating 7 ^ and t k is to use FFT (Fast 
Fourier Transform), which is both computationally efficient and robust to model 
errors. However, FFT is known for its high sidelobes and poor resolution. Many 
different types of windows may be applied to z to reduce the sidelobes. Yet using 
windows with FFT further reduces the already poor resolution of FFT. 

When the number of scattering centers K of a radar target is small, one may 
use many existing frequency estimation methods, such as the computationally and 
asymptotically (for high signal-to-noise ratio) statistically efficient MODE method 
[ 1 , 2, 3], to estimate the t k . Once t k are estimated, the straightforward least-squares 
fitting method may be used to estimate the corresponding 7 ^. 

However, there are two cases where the above parametric modeling approach may 
result in poor performance. First, for a complicated radar target such as an airplane, 
the number of scattering centers A may be very large and may even be larger than 
the number of samples M in z. For this case, the frequency estimation methods 
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cannot be used since to uniquely determine the time delays t k , even in the absence of 
the additive noise n 0 , the K must be less than (M - l)/2 [4, 5]. Second, when two 
or more scattering centers are very closely spaced, the frequency estimation methods 
may not be able to resolve all of the scattering centers. When this result occurs, the 
least-squares fitting method yields poor RCS estimates. As a result, the estimated 
range signature, i.e., 'Jk versus t k , may be distorted. 

We consider below a more robust parametric data model. To estimate the RCS 
7(f) of a scattering center whose range corresponds to a time delay f, we model the 
received data z as 

Z = 7(f)f(f) + n(f), (5-3) 

where f(f) is defined in ((5.2)) with t k replaced by t and n(f) denotes the additive 
noise. The n(f) is assumed to be circularly symmetric complex Gaussian random 
vector with zero-mean and unknown covariance matrix Q(i)- The unknown Q(f) 
models both the thermal noise and the interferences from other scattering centers. 
The Q(f) is assumed to be a positive definite and Toeplitz matrix. 

The problem of interest herein is to determine the range signature 7 (f) from z. 
We remark that to compute a discrete range signature with M s samples, we set 
f = 27r(m - 1 )/M s for the mth sample. We shall show below that the M s RCS 
estimates in the signature can be computed in parallel. 

5.2.2 The APES Algorithm 

We present below an algorithm that may be used to estimate the unknown 7(f) 
in model ((5.3)). We note that the amplitude and phase of 7 (f) may be considered as 
the amplitude and phase of the complex sinusoid with frequency f described by f(t). 
Thus we refer to our estimator below as the Amplitude and Phase Estimation of a 
Sinusoid in unknown colored noise , or simply the APES , method. 
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The APES algorithm is obtained by mimicking the computationally and asymp- 
totically statistically efficient large sample maximum likelihood estimator obtained m 
Appendix A for estimating the complex gain of a signal with known waveform and 
known steering vector in unknown colored noise. To obtain the APES algorithm, we 
first divide the measurement vector z into N overlapping subvectors z(n) of dimen- 
sion J\d x 1, where TV" = A/” — JVI T 1. The nth subvector z (n) contains the nth to 
the (n + M — 1 ) th element of z. Let n(f,n) be formed from n(t) in the same way 
z(n) are formed from z. The n(<, n) are assumed to be circularly symmetric complex 
Gaussian random vectors with zero-mean and the same unknown covariance matrix 
Q ,(*), which is a submatrix of Q(i). Then z(n), n = 1, 2, • • • , N, may be written as 

z(n) = 7 (i)a o{t)y{t, n) + n(f , n), (5.4) 

where ao(£) is referred to as the steering vector and has the form 

a 0 (f) — 1 , ( 5 - 5 ) 

and 

y(t,n) = e J ^ n-1)t , n = 1,2, • • • , N. (5.6) 

Let J be the M x M exchange matrix (with ones on the antidiagonal and zeros 
elsewhere). Let 

z(n) = Jz*(iV — n + 1), (5.7) 

where (•)* denotes the complex conjugate. Let 

n(t,n) = Jn*(t, N - n + 1), (5.8) 

where n(f,n) are assumed to have the same statistics as n(f,7i). Then z(n) may be 
written as 

z(n) = 7 '(t)/?(i)a 0 (<)y(t, n) + n (t, n), (5.9) 
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where 


(5.10) 


/?(() = e- i|w - 1) '. 

Let 

hy{t) = i (5.11) 

^ ’ 71 = 1 

where (•)^ denotes the complex conjugate transpose for matrices and complex con- 
jugate for scalars. Let r yy (t), Rzz, Rzz, and tzy{t) be defined similarly as ?zy(0- 
Note from ((5.6)) that r yy (t) = 1. Note also that as may be seen from ((5.11)), 
the mth element of h y (t) is the discrete-time Fourier transform of the mth sequence 
{z m (l), • • ■ , z m (N)} divided by N, where z m (n) denotes the mth element of z (n). 
Thus to compute a discrete range signature with M s samples, where M s is a power of 2, 
the h y (t), and similarly the may be computed with the FFT method and zero- 

padding. The amount of computations required for each sequence is 0[M S log a (Af,)]. 
Since the 2nd through the jVth elements of the mth sequence are the 1st through the 
(N - l)th elements of the (m + l)th sequence, the FFT of the (m + l)th sequence 
may be updated from the FFT of the mth sequence with 0(M S ) operations. 

Let b(i) = 7 (<)ao(<). Applying the results in Appendix A to z (n), we obtain the 
estimate of b(£) as 

£(0 = hy{t)- (5-12) 

Similarly, let b(<) = 7*(t)/?(f)a 0 (t)- Applying the results in Appendix A to z(n), we 
have 

B(<) = h y (t)- (5.13) 

By applying the results in Appendix A to both z (n) and z (n), the covariance matrix 
Qx(i) may be estimated as 

Ql(t) = ^ — ^zy(0^ly(0 "b R-zZ — Rzy(0^2y(0] 

= \ | Rzz + Rzz - [ b(f) 5(<)][b(0 6(0 ]*} 
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(5.14) 


= i[R + B(()B H (i)]. 

Then the estimate of 7 (t) may be computed as (see Appendix A) 

1 [ aTOQT^lM ± 

} 2 [ a^(t)Qr 1 (0ao(0 

Note that since only B(t) is a function of the range related time delay t, the Q^(t) 
can be computed more efficiently with 

Q-'(i) = 2 (r -1 - R-‘B(t) [6"(t)ft-'6(i) -l] _1 B"(()R-'} • (5.16) 

We note that computing R -1 requires 0(M 3 ) computations. Once R 1 is obtained, 
computing Q ^(t) for each t requires 0(M 2 ) computations. Computing 7 (t) for each 
t also requires 0(M 2 ) computations. 

Then the APES estimator for estimating q(t) may be summarized as follows: 

Step 1: Compute b (t), b (t), and Qi(t) with ((5.12)), ((5.13)), and ((5.16)), respec- 
tively. 

Step 2: Determine q(t) with ((5.15)). 

We remark that the vectors h(t,n) and h(t,n) are not independent of each other 
since they are formed as the overlapping subvectors of n(t). Yet since these vectors 
are not completely correlated with each other, it can be shown that as N goes to 
infinity, iz y {t), Rzz, Rzz, and r 2y (t) are all consistent estimates. As a result, the 
7 (t), for all possible t , is also a consistent estimate of 7(2). 

We also remark that the APES algorithm only requires simple matrix multiplica- 
tions and matrix inverses. Moreover, many computations of the APES algorithm can 
be done in parallel. It is thus possible to achieve real-time range signature estima- 
tion by implementing the APES algorithm with parallel processors and/or specially 
designed hardwares. 


(5.15) 
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We note that for the special case of M = 1, the 7 (f) in ((5.15)) can be shown to 
be the Fourier transform of z. Thus the Fourier transform method is a special case 
of the APES method. 

Finally, the parameter M has the following effects on the APES performance. 
Note first that the larger the M, the larger the dimension of Qi(t), and thus the 
better the modeling of the interferences. We shall show in the following subsection 
that the larger the Af, the better the resolution of the APES method. On the other 
hand, the larger the ibT, the smaller the N = M — M Thus increasing M increases 
the variance of Qi{t) since Rzzj Rzz> a nd *zy{t) are poorer estimates for larger 

M. We shall show in the following subsection that for very large M, the variance of 
7 (£) may increase. Also, increasing M increases the amount of computations needed 
by the APES method. 

5.2.3 Numerical and Experimental Results 

We present below both numerical and experimental examples showing the perfor- 
mance of the APES algorithm. In the simulated numerical examples, the root-mean- 
squared errors (RMSEs) of the APES estimator are obtained with 100 independent 
Monte-Carlo trials and are compared with the corresponding Cramer-Rao bounds 
(CRBs) derived in Appendix A. 

We first use simple numerical examples to illustrate the performance of the APES 
estimator. Consider first an example where the true range signature is shown in 
Figure 5 . 1 (a). We assume that there are 66 scattering centers in the signature. The 
number of data samples in the received data vector z in ((5.3)) is assumed to be 
M = 128. The additive thermal noise in z is assumed to be a zero-mean white 
Gaussian random process with variance 1 . Figure 5 . 1 (b) shows the estimated range 
signature (amplitude only) with the FFT method, which is also equivalent to the 
APES method with M = 1. Note that due to the large sidelobes of the FFT method, 
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the fourth scattering center is not discernible in Figure 5.1(b). Figure 5.1(c) shows the 
estimated range signature with the APES method with M = 48. Comparing Figures 
5.1(b) and (c), we note that using APES method with M = 48 yields a much better 
range signature estimate, which has much reduced sidelobes and better resolution for 
the large scattering centers, than using the FFT method. For the 62 very small and 
very closely spaced scattering centers between 0.97 and 2.47, the APES method with 
M = 48 tends to suppress them. Figure 5.1(d) shows the estimated range signature 
obtained when using the MODE with the least-squares fitting (MODE-LSF) method 
by assuming that there are 48 complex sinusoids plus white noise in the received data 
vector z. We note that since the data model is incorrect for the MODE-LSF method 
and some of the estimated scattering centers are very close to each other, using the 
MODE-LSF method may yield very large false peaks. For this case, they occur near 
t = 2.5. (Note that this type of large false peaks do not occur in every Monte-Carlo 
simulation. They occur in about 20% of the Monte-Carlo simulations.) The presence 
of these large false peaks makes the MODE-LSF method less preferable even than the 
FFT method. Figure 5.1(e) shows the estimated range signature obtained by using 
FFT with Kaiser window and shape parameter 4. Comparing Figures 5.1(b), (c), and 
(e), we note that using FFT with Kaiser window can reduce the FFT sidelobes, but 
the already poor resolution of the FFT method is made poorer due to the windowing. 

We now consider the effects of M on the performance of the APES estimator. 
Figure 5.2 shows the RMSEs of 7 (t = 0.2332) (the first scatterer in Figure 5.1(a)) 
obtained with the APES estimator as a function of M. We note from Figure 5.2 
that the APES estimator with a proper M can give much more accurate estimates of 
7(f = 0.2332) than the FFT method, which is equivalent to the APES method with 
M = 1 . With a proper M, the performance of the APES estimator can be close the 
corresponding CRB, which is also the best unbiased performance that can be achieved 
by an estimator. Note also that for very large M, the RMSE of 7 (i = 0.2332) 
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may increase due to the poor estimate of Qi(<)- To achieve the best estimate of 
~f(t = 0.2332), the M should be within 20 < M < 64. 

Figures 5.3(a) and (b) show the range signature estimates obtained with the APES 
method when M = 28 and M = 64, respectively. We note from Figure 5.2(a) that 
using the APES method with any M within 20 < M < 64 yields similar RMSEs for 

= 0.2332). Yet Figure 5.3 shows that the larger the M within 20 < M < 64, the 
better the resolution of the estimated range signature. The larger the M, however, 
the larger the amount of computations needed by the APES method. 

Consider next an example where the true range signature is shown m Figure 5.4(a). 
We assume that there are 75 scattering centers in the signature. The number of data 
samples in the received data vector z in ((5.3)) is assumed to be M = 128. The 
additive thermal noise in z is assumed to be a zero-mean white Gaussian random 
process with variance 10. Figure 5.4(b) shows the estimated range signature with the. 
FFT method. Note that due to the large sidelobes of the FFT method, the two small 
scattering centers that are to left and right of the first group of the large scattering 
centers are not discernible in Figure 5.4(b). Figure 5.4(c) shows the estimated range 
signature with the APES method with M = 48. We note that the two small scattering 
centers that are to left and right of the first group of the large scattering centers are 
resolved in Figure 5.4(c) since the sidelobes are very small for the APES method 
with M = 48. Comparing Figures 5.4(b) and (c), we note that neither the FFT nor 
the APES method with M = 48 can resolve the three groups of very closely spaced 
scattering centers. Both methods tend to combine the RCSs of very closely spaced 
scattering centers together to yield a large peak in the range signature since these 
scattering centers have the same phase (see Figure 5.4(a)). Figure 5.4(d) shows the 
estimated range signature with the MODE-LSF method by assuming that there are 48 
complex sinusoids plus white noise in the received data vector z. We note that since 
the data model is incorrect for the MODE-LSF method and some of the estimated 
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scatterer locations are very close to each other, using the MODE-LSF method may 
yield large false peaks. For this example, they occur near t = 1.7. Figure 5.4(e) shows 
the estimated range signature obtained by using FFT with Kaiser window and shape 
parameter 4. Comparing Figures 5.4(b) and (e), we note that using FFT with Kaiser 
window further widens the large peaks in 5.4(b). 

Finally, we apply below the APES method to the experimental data measured by 
a ground-to-air radar. (All other information about the radar and the radar target is 
not available for public release.) The measured data consists of M = 12S samples and 
is also degraded with the zero-mean white Gaussian noise with variance 0.1. Figure 
5.5(a) shows the target range signature of an aircraft obtained with the FFT method. 
Figure 5.5(b) shows the target range signature of the aircraft obtained by using the 
FFT method with the Kaiser window and shape parameter 4. Figure 5.5(c) shows the 
target range signature obtained with the APES method with M = 48. Comparing 
Figures 5.5(a) and (c), we note that using the APES method with M = 48 yields 
much lower sidelobes than using the FFT method. Comparing Figures 5.5(b) and 
(c), we note that using the APES method with M = 48 yields better resolution and 
sharper peaks than using the windowed FFT. Figure 5.5(d) shows the target range 
signature obtained when using the MODE-LSF method and assuming that there are 
48 complex sinusoids plus white noise in the measured data. We note that using the 
MODE-LSF method can yield poor target range signatures due to inaccurate data 
models. 

5.3 Using the APES Algorithm for SAR Imaging 

In this section, w r e extend the parametric model we proposed in the previous 
section to synthetic aperture radar (SAR) imaging. We also describe how the APES 
algorithm presented in the previous section may be used for SAR imaging. 
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5.3.1 Problem Formulation 

In synthetic aperture radar (SAR) or inverse SAR (ISAR) imaging [6, 7], the 
radar usually transmits linear frequency modulated (chirp) pulses. Upon receiving 
each pulse returned by an object being imaged, the radar demodulates the pulse by 
mixing the pulse with a reference chirp signal and low-pass filter the mixed signal. 
As a result, the scattering centers of the object at different ranges correspond to 
different frequencies of the output of this demodulation operation. Since either the 
radar or the object is moving or rotating, the pulses received at different angles 
between the radar and the object are used to form a synthetic aperture. After Polar- 
to-Cartesian interpolation, the scattering centers of the object at the same range but 
different cross-ranges correspond to different (Doppler) frequencies over the synthetic 
aperture. 

Let Z denote an M x x M 2 matrix containing the samples of the demodulated and 
Polar-to-Cartesian interpolated data, from which a SAR image is computed. A simple 
nonparametric method of computing the SAR image is to use the two-dimensional 
(2-D) FFT, which is again both computationally efficient and robust to model errors. 
However, FFT is known for its high sidelobes and poor resolution. Many different 
types of windows may be applied to Z to reduce the sidelobes. Yet using windows with 
FFT further reduces the already poor resolution of FFT. Thus to use 2-D windowed 
FFT to obtain high resolution SAR images, a large radar bandwidth and a large 
synthetic aperture are needed. One of the disadvantages of this requirement is that it 
makes it difficult to use synthetic aperture radar technology for wide-area surveillance. 

To achieve high resolution with limited radar bandwidth and synthetic aperture, 
many parametric and nonparametric spectral estimation methods have been used for 
SAR imaging [S, 9, 10, 11, 12, 13, 14]. For example, the Z may be written as the 
sum of 2-D complex sinusoids in additive noise. The 2-D frequencies of the sinusoids 
may be estimated with computationally and asymptotically (for high signal-to-noise 
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ratio) statistically efficient 2D-MODE method [14]. The complex amplitudes of the 
sinusoids may then be estimated with the straightforward least-squares fitting [9], 
However, similar to the case of the high resolution range estimation, this parametric 
modeling approach may result in poor performance when the number of sinusoids is 
too large to uniquely determine the parameters of the sinusoids [14] or when two or 
more sinusoids are so closely spaced that 2-D frequency estimation methods cannot 
resolve all of them. 

We consider below a more robust parametric data model, which is similar to the 
one proposed in the previous section. Let 7(^11 i 2 ) denote the complex intensity of the 
(z’,j)th pixel of an M sl x M s2 SAR image we intend to form, where U = 2 tt(z - l)/M sl 
and t 2 = 2t v(j - 1 )/M s2 . To estimate 7(^1, t 2 ), we model the received data Z as 

Z = 7(fi,i2)f(ii)®f T (fa) + N(f 1 ,*2), ( 5 - 17 ) 

where <S) denotes the Kronecker matrix product, f(^), i = 1>2, is defined in ((5.2)) 
with f/t replaced by t{ and M replaced by Mi , and N(ii,i 2 ) denotes the additive 
noise. The vec[N(ii, t 2 )]> where vec(-) denotes stacking all columns of a matrix into a 
single column vector, is assumed to be a circularly symmetric complex Gaussian ran- 
dom vector with zero-mean and unknown covariance matrix Q(ii>f 2 )- The unknown 
Q(*i,f 2 ) models both the thermal noise and the interferences from other pixels of 
the SAR image. The Q (ti,t 2 ) is assumed to be a positive definite and block-Toeplitz 
matrix whose blocks are also Toeplitz matrices. 

The problem of interest herein is to determine the range signature 7(^1 5^2) from 
Z. We remark that the complex intensities 7(^1, f 2 ) of a U pixels in a SAR image can 
be computed in parallel. 
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5.3.2 SAR Imaging with the APES Algorithm 

The APES algorithm may be applied to the rows (columns) and then to the 
columns (rows) of the data matrix Z to form the SAR image. We referred to this 
approach as the one-dimensional APES or the 1-D APES algorithm. We also present 
below a two-dimensional APES or the 2-D APES algorithm. We shall show in the 
next subsection that the 2-D APES algorithm may be computationally more intensi\e 
but may yield better SAR images than the 1-D APES algorithm. 

To obtain the 2-D APES algorithm, we first divide the data matrix Z into NiN 2 
overlapping submatrices Z (n) of dimension M\ x M 2 (in raster sequence), where 
JVi = Mi -M x A - 1 and N 2 = M 2 - M 2 + 1. Let N(t u t 2 ,n) be formed from N(h,t 2 ) 
in the same way Z(n) are formed from Z. Let 

h(ti,t 2 ,n) = vec [N^i,^,™)] . (5.18) 

The n(ti,t 2 ,n) are assumed to be circularly symmetric complex Gaussian random 
vectors with zero-mean and the same unknown covariance matrix Qi(ti,t 2 )) which 
has the same form as Q(ti,t 2 ). Let 

z(n) = vec [Z(n)j . (5.19) 

Then z(n), n = 1,2, ■ ■ • , NiN 2 , may be written as 

z(n) = 7 (ti,t 2 ) a o(^i^ 2 )j/(^i)^ 2 ,”) + n(fi,f 2 ?*i)j (5.20) 

where a. 0 (ti,t 2 ) is referred to as the steering vector and has the form 

ao(^i) I 2 ) — a i(^i) ® a 2 (^ 2)1 (5.21) 


where 


i(U) = 


Jt < 




-1 i 

■ 2 = 12 

5 1 


(5.22) 
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and y(£i,£ 2 ,n) is the nth element of the vector yo^i,^)* 

yofaifc) = yi(*i) ®y2(<2), 


(5.23) 


with 

y,(£,)=[l e J( ‘ ] , * = 1,2. (5.24) 

Let J be the M\M 2 x M\M 2 exchange matrix. Let 

z(n) = Jz*(iV x fV 2 - n + 1), (5.25) 

Let 

n(t\, £ 2 , n) = Jn*(£ x , £ 2 , NiN 2 — n + 1), (5.26) 

where n(£ 1 ,£ 2 ,n) are assumed to have the same statistics as n(£i,£ 2 ,n). Then z(n) 
may be written as 

z(n) = 7 *(£i, £ 2 )/3(£i, £ 2 )ao(£i, £ 2 )j/(£i, £ 2 , n) + n(£i,£ 2 ,n), (5.27) 


where 

P(t u t 2 ) = (5.28) 

With these new notation definitions, the remaining steps of computing 7 (£i,£ 2 ) in 
the 2-D APES algorithm are similar to those of computing -y(£) in the APES algorithm 
presented in the previous section. Let 


h y {tuh) = 


1 N x N 2 
£Viiv 2 n= i 


(5.29) 


and r ay (<i,* 2 ) be defined similarly. Then similar to using FFT to compute f iy (t) and 
fg 9 (£) in the previous section, the 2-D FFT method can be used to compute fzy(£ij £ 2 ) 
and f , 2 !/ (£i,£ 2 ) more efficiently. 

We remark that for the special case of M\ = M 2 — 1» the 7(£i,£ 2 ) obtained 
with either the 1-D or the 2-D APES algorithm can be shown to be the 2-D Fourier 
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transform of Z. Thus the Fourier transform method is again a special case of the 
APES methods. 

We also remark that both the 1-D and 2-D APES algorithms only require simple 
matrix multiplications and matrix inverses. Moreover, many computations of the 
APES algorithms can be done in parallel. It is thus possible to achieve real-time 
SAR image formation by implementing the APES algorithm with parallel processors 
and/or specially designed hardwares. 

5.3.3 XPATCH Examples 

We present below simulated examples showing the performance of the 1-D and 
2-D APES algorithms for SAR imaging. The data we use are generated by XPATCH, 
which is a computational electromagnetics software [15]. 

Figure 5.6(a) shows the SAR image (the logarithm of the amplitude) of a firetruck 
obtained with the 2-D FFT method when Mi = M 2 — 128. The firetruck is simulated 
with the XPATCH software [15] and is added with the zero-mean white Gaussian 
noise with variance 0.1. Figure 5.6(b) shows the SAR image obtained by using the 
2-D FFT method with the circularly symmetric Kaiser window and shape parameter 
6. Figure 5.6(c) shows the SAR image obtained by using the 1-D APES algorithm 
with M = 48. Comparing Figures 5.6(a), (b), and (c), we note that using the 1- 
D APES algorithm with M - 48 yields lower sidelobes and better resolution than 
the 2-D FFT methods. Figure 5.6(d) shows the SAR image obtained by using the 
2D-MODE algorithm with the least-squares fitting (2D-MODE-LSF) and assuming 
that there are 121 two-dimensional complex sinusoids in additive white noise. We 
note that the 2D-MODE-LSF algorithm performs poorly and many features that are 
important for target identification, such as the target pixels to the upper right of the 
brightest pixels in Figure 5.6(b), are left out in Figure 5.6(d) by the 2D-MODE-LSF 
algorithm. 


Figure 5.7 is similar to Figure 5.6 except that only 1/2 of the bandwidth and 
1/2 of the synthetic aperture are used for SAR imaging. Further, Figure 5.7(c) is 
obtained with the 1-D APES algorithm with M = 24. Figure 5.8 is also similar to 
Figure 5.6 except that only 1/3 of the bandwidth and 1/3 of the synthetic aperture 
are used for SAR imaging. Further, Figure 5.8(c) is obtained with the 1-D APES 
algorithm with M = 16. Figure 5.8(d) is obtained with the 2-D APES algorithm with 

= M 2 — 16. Figure 5.8(e) is obtained similarly as Figure 5.7(d). We again note 
that using the APES algorithms yield lower sidelobes and better resolution than the 
2-D FFT methods. The APES algorithm is also more robust than the 2D-MODE-LSF 
method. Finally, we note from Figures 5.8(c) and (d) that the 2-D APES algorithm 
gives better performance than the 1-D APES algorithm at the cost of much more 
computations. 

5.4 Conclusions 

We have presented a robust parametric data model for estimating high resolution 
range signatures of radar targets and for forming high resolution synthetic aperture 
radar (SAR) images. We have described how the APES (Amplitude and Phase Es- 
timation of a Sinusoid in unknown colored noise) algorithm can be used to estimate 
the range signatures and to form SAR images. With both numerical and experi- 
mental examples, we have shown that our modeling and estimation approach yields 
better resolution and lower sidelobes than the conventional nonparametric FFT (fast 
Fourier transform) method. We have also shown that our approach is more robust 
than modeling the radar data as a certain number of complex sinusoids in noise and 
estimating the frequencies, amplitudes, and phases of the sinusoids with one of the 
best sinusoidal parameter estimation methods. 
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Appendix A — Large Sample Maximum Likelihood Estimator and Cramer- 
Rao Bound 

We consider below the estimation of the complex gain of a signal with known 
waveform and known steering vector in unknown colored Gaussian noise. Let x(n), 
n = 1,2, • • • , N, denote some received data vectors of dimensions M x 1. Let y(n), 
n — l ; 2, , N, denote a known waveform. Assume that the received data vectors 

x(n) can be written as 

x(n) = 7 a 0 y(n) + n(n), (5.30) 

where n is the additive noise vector of dimensions M x 1 and ao is referred to as the 
steering vector and is known. The noise vectors n(n) are assumed to be circularly 
symmetric complex Gaussian random vectors with zero-mean and arbitrary covariance 
matrix Q and are temporally white, i.e., 

E[n(i)n ff (j)] = QS itj , (5.31) 

where (-) H denotes the complex conjugate transpose and is the Kronecker delta. 

We consider below a large sample (JV > 1) maximum likelihood (ML) estimator 
for estimating 7 from x(n), n = 1,2,---,AL It is easy to show that an exact ML 
estimator requires a multidimensional search over the parameter space and is com- 
putationally burdensome. We shall present below a large sample ML estimator that 
is both computationally and asymptotically (for large N ) statistically efficient. The 
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approach we use to derive the large sample ML estimator is similar to the one in 
[16, 17], which is devised for a different problem. 

The log-likelihood function of the received vectors x(n), n — 1, 2, • • • , N, is pro- 
portional to (within an additive constant): 

- In |Q| - tr |Q _1 ^ Et x ( n ) - by(n)][x(n) - by(n)] H | , 
where ] • | denotes the determinant of a matrix and 

b = 7ao- (5.32) 


Consider first the estimate of Q and the unstructured estimate of b. It is easy to 
show that 

Q = T7 £[x(n) - by(n)][x(n) - by{n)] H , (5.33) 

^ 71=1 

and b may be obtained by minimizing the following cost function 


F = 


T7^[ x (”)- b y( n )][ x ( n )- b 2/( n )] H ' 

71 = 1 


(5.34) 


Let 


and 


? 


yx — 


r yy “ 


T7 E 

Iy 71 = 1 

T7 J2y( n )y H ( n )- 

iV 71=1 


Let R xx be defined similarly as r yy . Then let [18] 


(5.35) 

(5.36) 


G = 


Et x ( n ) ~ b y( n )H x ( n ) - b ^( n )] ff 

^ 71=1 

R** - bf s x - + bf„b" 


u L yx l yy 


yy 


~H --1 
r yx r yy 


H 


+ Rxx - r yx r y ^r yx . 


(5.37) 

(5.38) 

(5.39) 
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Since fyy is a positive number and the second and third terms in ((5.39)) do not 
depend on b, it follows that 


G > G 


b=b> 


(5.40) 


where 


u -H — 1 

b = Vyy • 


(5.41) 


Since the whole sample covariance matrix G is minimized, the unstructured estimate 
b of b in ((5.41)) will minimize any nondecreasing function of G including the deter- 
minant of G, which is-jF in ((5.34)). It is easy to see that b is a consistent estimate 
of b. 

By using ((5.41)) with ((5.33)), the Q may be rewritten as 


Q — R-XX — TyxTyy fyx- 


(5.42) 


It is easy to see that Q is a consistent estimate of Q. 

Let us now consider the structure of b. The cost function in ((5.34)) may be 

rewritten as 


F = 


Rxx - bf yx - ff x b ff + br. 




yy ‘ 


R 

Q 


XX - br yy b H + (b - b)r yy (b - b) H 

I |l + Q _1 (b — b)r yj ,(b — b) H I . 


(5.43) 

(5.44) 

(5.45) 


The ML estimate of 7 may be obtained by minimizing F in ((5.45)) or equivalently 
In F. In [16, 17], we have shown that minimizing InF, with F defined in ((5.45)), is 
asymptotically (for large N ) equivalent to minimizing 


F 1 =tr[r„„(b-b)' , Q-’(b-b)]. (5.46) 


Using b = qao and minimizing in ((5.46)) with respect to 7 yields 


7 


agQ^b 

a o Q _1 ao 


(5.47) 


106 



The large sample ML estimator for estimating 7 may be summarized as follows: 


Step 1: Compute b and Q with ((5.41)) and ((5.42)), respectively. 


Step 2: Determine 7 with ((5.47)). 

We remark that since 7 is a consistent and large sample realization of the ML esti- 
mate, it follows that 7 is asymptotically (for large N) statistically efficient according 
to the general properties of ML estimators [19]. Using the results in [16, 17], we can 
show that the asymptotic (for large N ) distribution of 7 is complex Gaussian with 
mean 7 and variance 


Var[ 7 ] = 


1 

EJLi Q _1 aoy(») ’ 


(5.48) 


which is also the CRB. 
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Figure 5.1: Range signature estimates (solid lines in (b) - (e)) compared to the true 
range signature (dashed lines in (b) — (e)) when the white noise variance is 1 and 
M = 128. (a) True range signature, (b) FFT or APES with M = 1. (c) APES with 
M = 4S. (d) MODE-LSF by assuming 48 complex sinusoids, (e) FFT with Kaiser 
window and shape parameter 4. 
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Figure 5.4: Range signature estimates (solid lines in (b) - (e)) compared to the true 
range signature (dashed lines in (b) - (e)) when the white noise variance is 10 and 
M — 128. (a) True range signature, (b) FFT or APES with M = 1. (c) APES with 
M = 48. (d) MODE-LSF by assuming 48 complex sinusoids, (e) FFT with Kaiser 
window and shape parameter 4. 
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Figure 5.6: S A R images of an XPATC'H simulated firetruck when the additive white 
noise variance is 0.1 and d/j = d/j = 12S. (a) 2-D FFT. (b) 2-D FFT with circularly 
symmetric Kaiser window and shape parameter 6. (c) 1-D APES with M — -IS. (d) 
2D-MODF.-LSF by assuming 121 complex sinusoids. 
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Figure 5.7: SAR images of an XPATC'H simulated nretruck when the additive white 
noise variance is 0.1 and - M 2 = 128. Further. 1/2 of the bandwidth and 1/2 
of the synthetic aperture of that in Figure 5.6 are used, (a) 2-D FFT. (b) 2-D FFT 
with circularly symmetric Kaiser window and shape parameter 6. < c 1-D APES with 
-1/ = 2'. id' 2D-MODE-LSF by assuming 121 complex sinusoids. 
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Figure 5.3: SAR images ot an XPATCH simulated firerruck when the additive white 
noise variance is 0.1 anci M x = M : = 128. Further. 1/3 of the bandwidth and 1/3 
ot the synthetic aoerture of that in Figure 5.6 are used, (a) 2-D FFT. (bi °-D FFT 
with circularly symmetric Kaiser window and shape parameter 6. < c i !-D APES with 
3/ = ib. f.l: 2-D APES with V, = .U 2 = 16. (e) 2D-M0DE-LSF by assuming 121 
c o m p leu s i n 1 iso i ri > . 
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Figure 8: (e) 
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